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of a thermoacoustic prime mover with periodic boundary conditions. There are five 
significant aspects to this research: (1) using DeltaE to analyze an annular prime mover, 
(2) developing an entirely new analysis program using MATLAB, (3) designing, building, 
and experimentally investigating a single stack, annular prime mover, (4) experimentally 
investigating a constricted, single stack annular prime mover, (5) predicting the 
performance of a two stack annular prime mover. The major conclusions are: (1) A single 
stack annular prime mover will not reach onset because the eigenmodes of the system do 
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because the constriction forces dominating boundary conditions that alter the eigenmodes. 
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I. INTRODUCTION 


The subject of this dissertation is the investigation of thermoacoustic prime movers in 
an annular geometry. What sets this research apart from previous work in prime movers is 
the nature of the boundary conditions. To our knowledge this work is the first thorough 
investigation of thermoacoustic prime movers with periodic boundary conditions. The 
primary conclusion drawn from this work is that a full understanding of the eigenmodes of 


the system are required to design a functioning annular prime mover. 


Thermoacoustic heat transport is a process through which an acoustic field generates, 
or is generated from, a flow of heat. Thermoacoustic engines are of two types: heat pumps 
and prime movers. A thermoacoustic heat pump utilizes a standing wave to transport heat 
along the boundary of a plate situated in the standing wave. The acoustically generated heat 
flow produces a temperature gradient across the stack. In other words, acoustic energy is 
converted into stored thermal energy. In contrast, a thermoacoustic prime mover produces 
acoustic work, by accepting heat from a high-temperature source and transferring it to a 
low-temperature sink. In this dissertation a prime mover in an annular resonator is 


investigated both theoretically and experimentally. 


To ulustrate how a prime mover operates, two conventional prime mover 
configurations are described. One type (as is shown in Fig. 1.1) is a rigid-rigid standing 
wave tube which contains a stack of plates called the prime mover stack (or, simply, the 
stack), which is in thermal contact with two heat exchangers. In operation, a temperature 
difference is imposed between the two heat exchangers, resulting in a temperature gradient 
along the stack. When the temperature gradient in the stack 1s sufficiently large, the gas in 


the resonator oscillates spontaneously at a certain frequency, with pressure antinodes at the 


closed ends. The inception of spontaneous oscillation is known as onset. Once onset is 
reached, the acoustic amplitude in the resonator grows rapidly, typically reaching several 


percent of the mean gas pressure. 


Rigid Boundary Ambient Heat Exchanger Hot Heat Exchanger 





Ambient End Prime Mover Stack Hot End 


Figure 1.1 A typical rigid-rigid prime mover configuration. 


The second prime mover configuration, called the Hofler tube, is a mgid-open 
resonator (as is shown in Fig. 1.2). In 1983, Hofler designed and built this simple 
thermoacoustic oscillator to demonstrate some thermoacoustic phenomena to his doctoral 
committee at the University of California, San Diego. The sound radiated from the open 
end of the pipe is impressively loud, about 100 dB re 20 uPa a meter away (Wheatley et 
al., 1985; Swift, 1988). 


Ambient Heat Hot Heat 


Exchanger Exchanger 
! == = = 
Open | = —— Rigid 
Boundary , = === Boundary 
=: —<<— 
Ambient End Prime Mover Stack Hot End 


Figure 1.2 The Hofler tube, an example of a rigid-open prime mover. 


The condition necessary for a prime mover to reach onset is that the amount of stored 
thermal energy converted to acoustic energy must exceed the total amount of acoustic 
energy dissipated by thermal-viscous losses in the prime mover. This ability to reach onset 
is determined by the geometry of the prime mover, in particular, the position of the stack 
relative to the nearest pressure antinode. The hot end of the stack 1s closer to a pressure 
antinode than is the ambient end. In the two examples of conventional prime movers 
described above, the node and antinode positions are fixed by incorporating well-defined 
acoustic boundary conditions into the system. This is an important common feature for all 
conventional prime movers. The “built-in” boundary conditions also serve as convenient 
starting points for computational analysis of the performance of a prime mover. Because 
the relative positions of the dominating boundaries and the stack are fixed, the stack 
position is optimized for only one acoustic mode. Also, the optimal stack position for 
onset is not necessarily the optimal position for high amplitude performance. One of the 
initial motivations for studying annular prime movers was to see if it would self-optimize 


the stack location relative to the acoustic field. 


Unlike a conventional prime mover, an annular prime mover (as is shown in Fig. 1.3) 
does not have the typical dominating boundary condition to force a pressure or velocity 
node at any particular position. Design and analysis of a functional annular prime mover 


are thus made more difficult. 





Hot Heat Exchanger Ambient Heat Exchanger 


Hot End Ambient End 


Prime Mover Stack 


Figure 1.3 An annular prime mover configuration. 


In a uniform cross-section annular resonator, the fundamental longitudinal mode 
corresponds to the circumference being approximately equal to one acoustic wavelength. 
There are two degenerate orthogonal modes that satisfy this condition. However, these 
modes have no preferred orientation. The presence of a nonuniformity (for example, the 
stack) breaks the degeneracy resulting in a frequency-splitting. The low frequency mode 
will have a pressure node at the stack, and the higher frequency mode will have a velocity 


node at the stack. (In what follows, the word “frequency” will be dropped and the modes 


referred to, simply as “high” and “low”.) Neither of these modes _ supports 


thermoacoustics. 


Although in retrospect the answer now seems obvious, prior to a thorough study of 
the problem one might have argued that, once it becomes very strong, the thermoacoustic 
effect may dominate the situation and shift the orientation of the acoustic field to some 
preferable position. This suggests the possibility of the ana prime mover achieving 
onset. Even if the uniform cross section, single stack prime mover does not reach onset, 
other nonuniformities in cross section will have some effect on the spatial distribution of 
the acoustic field. Hence, by placing another constriction in the annulus, the previously 
mentioned pressure and velocity nodes may be displaced from the stack, providing the 
possibility of the annular prime mover reaching onset. These are some of the concepts that 


initiated this work. 


It should be noted that a rigid-rigid prime mover can be considered a limiting case of a 
constricted, annular prime mover. Al] the previous work on thermoacoustic prime movers 
have dealt with engines that have some sort of well-defined boundary conditions imposed 
upon them. There has been no detailed analysis of prime movers with periodic boundary 


conditions. 


A. BACKGROUND 


Although thermoacoustic phenomena have been observed for a long time, significant 
advances in practical thermoacoustics are relatively recent. The earliest example of a 


thermoacoustic prime mover is the Sondhauss tube (Sondhauss, 1850). Over 100 years 


ago, it was observed by glassblowers that a hot glass bulb attached to a cool glass tube 
sometimes emitted sound at the tip of the tube. Sondhauss quantitatively investigated the 
relationship between the pitch of the sound emitted and the dimensions of the tube. Lord 
Rayleigh explained the Sondhauss tube quantitatively in 1896, but no complete theoretical 
analysis of these phenomena was made before publication of the series of papers by 
Nikolaus Rott (Rott, 1969, 1980, 1983). Later, Wheatley (1985) and others at Los 


Alamos began the development of practical thermoacoustics devices. 


B. PAST RELEVANT WORKS 


Past works relevant to this research fall under three headings: (1) the thermoacoustic 


prime mover, (2) the acoustic Stirling Engine, and (3) the annular resonator. 
1. The Thermoacoustic Prime Mover 


The potential application of thermoacoustics as a heat-driven sound source has 
motivated theoretical developments (Rott, 1969; Wheatley and Cox, 1988; Swift, 1988). 
Work on prime movers has been focused on either rigid-rigid or rigid-open prime movers. 
Some examples are outlined below. Miglori and Swift (1988) constructed a thermoacoustic 
prime mover that used liquid sodium as its working fluid. A prime mover has also been 
used as a heat-driven sonar projector (Gabrielson, 1991). Swift (1992) built and analyzed 
a 5-inch thermoacoustic prime mover which, at its most powerful operating point, using 
13.8-bar helium, delivered 630 W to an external acoustic load. With minor modifications, 
this device was used by Swift (1994) to research the application of similitude to nonlinear 


thermoacoustics. 


The Naval Postgraduate School has been in the forefront of much thermoacoustic 
research. Atchley and the author first built and tested a heat driven prime mover (Lin and 
Atchley 1989), which generated a sound, at a temperature difference of 453 °C, with a 
peak acoustic amplitude of 7.9% of atmospheric pressure. Subsequently, other work was 
done on this device (Atchley 1992, 1993; Atchley and Kuo, 1994). In 1993, a 
thermoacoustic prime mover was constructed by Castro and Hofler to investigate the 
performance of heat exchangers and produced peak pressures of up to 20% of the mean 
pressure (Castro and Hofler, 1993). 

DeltaE, a program developed at Los Alamos National Laboratory by Ward and Swift 
(1994, 1996) has been applied to several cases in the design and analysis of thermoacoustic 
prime movers. DeltaE stands for Design Environment for Low Amplitude Thermoacoustic 
Engines. It solves Rott’s wave equation for a geometry defined by the user. For example, 
Swift (1992, 1996) used DeltaE to analyze the performance of a 5-inch thermoacoustic 
engine. Nessler (1994) used DeltaE to compare the performance of the pin stack with a 
conventional stack in a thermoacoustic prime mover (Swift and Keolian 1993). Yang 
(1995) and Meng (1996) also used DeltaE to predict onset of their thermoacoustic prime 
movers. More recently, DeltaE was applied to analyze the efficiency of a thermoacoustic 


prime mover with a pin stack (Gibson, Nessler, and Keolian, 1997). 
2. The Acoustic Stirling Engine 


Ceperley (1979, 1982 and 1985) has discussed acoustic engines using traveling 
waves. As recognized by Ceperley and several other authors (Swift, 1988, 1995; Hofler 
1988; Organ 1992), the phasing between pressure oscillations and velocity within a Stirling 
engine regenerator is the same as those of an acoustic traveling wave. Ceperley proposed 


the replacement of the function of the pistons in a Stirling engine by acoustic processes, to 


form what he called a traveling-wave heat engine. One of his concepts, a traveling-wave 
heat-driven refrigerator, is shown in Fig. 1.4. In this device, one regenerator and heat 
exchanger set functions as the prime mover, adding acoustic power onto the traveling wave 
as heat flows from a high-temperature heat source to a room temperature heat sink. The 
other set functions as a heat pump, using acoustic power from the traveling wave to 
transport heat from a low temperature to room temperature. His work was partly 
responsible for motivating the author to embark on this project. However, his analysis was 


overly simplistic and he never built a working device. 
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Figure 1.4 Ceperley’s design for a traveling-wave heat-driven refrigerator. The 


arrows show the direction of wave propagation. 
3. The Annular Resonator 


Previous work on annular resonators is also relevant to the study of the annular prime 


mover. The eigenmodes for electromagnetic wave propagation in an annular cavity are of 


fundamental interest in the radio-frequency (RF) heating of tokamak plasmas, besides 
having other applications in the microwave circuit theory (Cap and Deutsch 1978, 1980; 
Janaki and Dasgupta, 1990; Wu, 1992). 

Another category of relevant work is the frequency splitting phenomenon found in 
different fields. As early as 1969, Rudnick et al. made observations of a superfluid helium 
persistent current using the Doppler-shifted splitting of an: azimuthal resonant fourth-sound 
mode of the cylindrical resonator containing liquid helium (Rudnick et al. 1969). 
Heiserman investigated the persistent currents in superleaks in contact with bulk superfluid 
helium using the Doppler shifts of the acoustic modes of an annular resonator partially 
packed with a superleak, and with a simple gyroscopic technique (Heiserman, 1975). In 
plasma physics, in the presence of a finite plasma current, the axisymmetric 
magnetoacoustic wave resonance exhibits a frequency splitting for a finite annular mode 
number between oppositely directed traveling waves (Borg and Wit, 1991). A deformation 
of the cross section of a pinch of the plasma cross section also induces frequency shift and 
Splitting of circular modes (Ring, 1988). In the microwave and RF field, frequency 
splitting can also be found in a nonuniform microstrip ring resonator (Wested and 
Andersen, 1972), in an asymmetrically coupled microstrip ring resonators (Al-Charchafchi 
and Dawson, 1990) and 1n a microstrip rhombic resonator with a step discontinuity in one 


of the arms (Al-Charchafchi and Boulkos, 1990; Al-Charchafchi and Schreck, 1994). 


C. OUTLINE OF THIS WORK 


This work constitutes the first detailed theoretical and experimental investigation of a 
thermoacoustic prime mover with periodic boundary conditions. There are five significant 


aspects to this research: (1) using DeltaE to analyze an annular prime mover, (2) 


developing an entirely new analysis program using MATLAB, (3) designing, building, and 
experimentally investigating a single stack, annular prime mover, (4) experimentally 
investigating a constricted, single stack annular prime mover, (5) predicting the 
performance of a two stack annular prime mover. 

The major conclusions are: (1) A single stack annular prime mover will not reach 
onset because the eigenmodes of the system do not support thermoacoustic’growth. (2) A 
constricted annular prime mover will reach onset because the constriction forces dominating 
boundary conditions that alter the eigenmodes. (3) A two stack prime mover is predicted to 
reach onset because one of the eigenmodes of the symmetric system does support 


thermoacoustic growth. 
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Il. THEORY 


This chapter is divided into four parts. The first part is a brief review of the aspects of 
thermoacoustics pertinent to the annular prime mover problem. The review draws heavily 
from Swift’s work (Swift, 1988 and 1997), and is intended to introduce the concepts 
required to understand the techniques used to analyze the annular prime’ mover. One of the 
central questions to be addressed in this dissertation is the ability of an annular prime mover 


to reach onset. The dependence of quality factor Q of the prime mover on the temperature 


difference AT applied to the stack is a key measure of this ability. In this research, the 


predicted Q will be determined from the complex eigenfrequency of the prime mover. 
Therefore, following the review of thermoacoustics, the complex eigenfrequency will be 


introduced. 


Another important aspect of this dissertation is the numerical analysis of the annular 
prime mover. The techniques used will be discussed in the last part of this chapter. Some 
basic properties of the two-point boundary-value problem are discussed, followed by a 
discussion of the advantages and disadvantages of using DeltaE with the annular prime 
mover geometry. A MATLAB numerical analysis program is then described that is more 
tailored to the annular prime mover than is DeltaE. The validation of this program is 
discussed and the application of the MATLAB program to our annular prime mover is 


presented. 


This chapter concludes with a discussion of the end corrections for constrictions in the 


cross section of annular resonators. 


A. REVIEW OF THERMOACOUSTICS 


Thermoacoustics is essentially the study of the acoustics of a fluid in which there 
exists a temperature gradient. The acoustics of such a fluid is contained in three differential 
equations which relate the complex acoustic pressure and volume velocity, and the 
temperature and enthalpy of the fluid. This method was first described by Wheatley et al. 
(Wheatley et al., 1983). The derivations of these relationships are summarized in this 
section. Once these concepts have been introduced, a description of the basic scheme used 
to analyze the performance of an annular prime mover is given. 

This section is essentially a summary of pertinent aspects of thermoacoustics given in 
the articles by G. W. Swift (Swift, 1988 and 1997). Swift’s quantitative thermoacoustic 
analysis is based on Rott’s wave equation and energy-flux equation, which result from the 
momentum, mass continuity, and energy equations in the acoustic approximation. The 
results are valid for arbitrary phasing between the acoustic pressure and velocity; in other 
words, for both standing-wave and traveling-wave systems. The results presented in this 
section serve as the basis for the numerical analysis used in this dissertation. 

First, the notation that will be used throughout this discussion is established. The 
acoustic wave is considered to propagate in the x direction within a duct of constant cross 
section. The usual complex notation is adopted for time-oscillatory acoustic quantities. 
Also, expansions to first order in the acoustic variables is assumed to suffice. Thus, the 


pressure, velocity and temperature can be written 


P(x)=p +P, el, (235) 


iGo eee y,z)ye. 5 (2.2) 


and 


2 


Cain) = 1 ee Taceyezne™ , (2.3) 


where i =-/—1 and u, is the x-component of y. The mean values (subscript m) are real, 
but the acoustic variables (subscript 1) are, in general complex, to account for the relative 
phasing of those oscillating quantities. It 1s assumed that the mean fluid velocity u,, = O. 
All the time dependence appears in the e“ term, with @ = 2 7 f being the angular 
frequency. For clarity, we consider here only the case of large solid heat capacity, so that 
the temperature of the solid material in the stack 1s simply T,,(x), independent of t, y, and z. 


The conservation of momentum of an incompressible viscous fluid is ( Landau, 1975) 


+ (v egrad)y = - <; grad p ys V*y, (2.4) 


|S 


where v is the velocity, p is the density, p is the pressure and p is the dynamic viscosity of 
the fluid in the resonator. This is called the Navier-Stokes equation. To develop a 
quantitative understanding including viscous effects, we begin by finding the y and z 
dependence of the velocity. Because of the viscous boundary layer in the y and z direction 
all other viscous derivatives can be neglected compared to pd*u,/dy* and p07u,/dz’. 


Equation (2.4) then reduces to 








' d O*u O°u 
iM p,u, = sim + I Se fi |. (2.5) 


Equation (2.5) is an ordinary differential equation for u,(y, z). The boundary 
condition at the solid surface is u, = 0. With this boundary condition Eq. (2.5) can be 


solved to yield 
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u(x, y, Z)= [1-h0, yh (2.6) 


m 


where /1,(y, z) depends on the specific geometry under consideration. The volume velocity 
is the spatial average over y and z of Eq. (2.6) and is thus 
TT A as d, 
U(x) = i (1 f,) 2, (2.7) 
OD) dx 


where f, is the spatial average of h,(y, z). For the case of parallel-plate geometry of the 


stack used in this work, it can be shown that (Swift 1988; Arnott et al. 1991) 


_ cosh[(1+2)y/d, ] 
MY) = cosh + Dy, 15,1" oe 
and 
tanh[(/ +i)y,/0, ] 
ee eee 2.9 
h id +2) y_/0,.] S” 


where 6, = J2H/P,@ is the fluid’s viscous penetration depth and y, is plate half-gap. 
Calculation of the oscillating fluid temperature is accomplished using only the first 
order terms in the general equation of heat transfer (Landau, 1975). The entropy is 
expressed in terms of the acoustic pressure and temperature of the fluid and the resulting 
ordinary differential equation is solved subject to the appropriate boundary conditions. In 


the presence of a temperature gradient along the duct, the solution is (Swift 1988 and 1997) 








; aT h 26h = 
18 (6, at [DN], 2 a i ll a (0) 
Hop Daj aye ee: = eg 


where O=c, [/k is the Prandtl number, and c, is the isobaric heat capacity per unit mass. 


The spatial average in (y, z)-direction of Eq. (2.10) can be written 
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where f, is the spatial average of h,. The functions h, and f, are the same as h, and f,, only 

with 6, replaced by 6, =./2K/p,,c,@ , the thermal penetration depth, where Kis the thermal 

conductivity of the fluid. The function f, and f, are very important.in thermoacoustics and 

describe how viscous and thermal processes affect the oscillatory velocity and temperature 
in the acoustic field. A graph of f, as a function of y/6, is shown in Fig. 2.1. 

The equation of continuity is (Landau, 1975) 


P+ Ve (pi) =0. (2.12) 


To first order, the x-component of Eq. (2.12) can be reduced to the form 


d(p,, %) 


1Mp, + = 


= 0. (28) 


Equation (2.13) can be rewritten in terms of volume velocity 





iMp, + Fi =a =a (2.14) 


gas 


Combining Egs. (2.10) and (2.13) yields 


tg | A 
= = SLM les + 7 


_— ef.) 1 Ta 
dx mc 1-f,)\1-o)T,, dx 


U, (2.1) 





where c is the speed of sound in the gas, and y is the ratio of isobaric to isochoric specific 


heats. 
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We pause to discuss the results derived in the last several paragraphs. Equation (2.7) 
shows that the gas which is much farther than 6, from the nearest solid surface experiences 
essentially no viscous shear and moves with a velocity that depends only on x. The gas 
that is much closer than 6, is nearly at rest. Gas that is between these extremes moves with 
a reduced velocity amplitude and significant phase change. The terms in Eq. (2.11) 
indicate that there are two contributions to the oscillatory temperature. The-first term is 
simply due to the adiabatic compressions and expansions of the fluid. The second term 
comes from the convection of gas parcels along a temperature gradient. The net 
temperature oscillation is just a linear combination of, these two effects. It should be 
pointed out that f. has the same functional dependence as y/6, as is portrayed in Fig. 2.1 
Co) ae ae 

Equation (2.15) has an easy physical interpretation. It indicates that dU, /dx comes 
from a density oscillation in the fluid. This density change can be caused by a pressure 


change or by convection of gas along the temperature gradient. 
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Figure 2.1 The real and imaginary parts and magnitude of f, as functions of the ratio 


y/6, for the case of parallel plate geometry. 


Next, following the derivation by Swift (1988 and 1997), the time averaged energy 


flux along the stack can be found as 


- We ce 
= 5 Re Ui ( al 
dT, 


ge a + of, )—*- (A K+A_.K “4 
2@ ree (_- ao”) IL - f, ; i "! dx gas soild ™ solid 


(2.16) 
aT, 


dx’ 





where the superscripts * represent the complex conjugate of a complex quantity, and K,,);, 
is the thermal conductivity of the solid. On the basis of Swift’s assumptions, in steady 
state for an engine without lateral heat flows to the surroundings, H, along x must be a 


constant throughout the stack. 


To numerically analyze the performance of thermoacoustic engines, it is recognized 
that Eqs. (2.7), (2.15) and (2.16) comprise a set of three second-order coupled complex 
differential equations in the five variables Re(p,), Im(p,), Re(U), Im(U;) and T_. These 
equations are the basis of both DeltaE and the MATLAB program presented later in this 
chapter. 

In order to obtain a more useful form of the differential equation for p,(x), Eq. (2.7) 


can be written as 


UIE UO? 


= A_(1- fi.) Ome.) (Qe 


In a duct region, the integration is governed by Eqs. (2.15) and (2.17) and the temperature 
gradient dT /dx is determined by the measured temperature profile. It is assumed that there 
is no temperature gradient across the heat exchanger. As a result, in a heat exchanger, 
integration is carried out by Eq. (2.17) and a simplification of Eq. (2.15) 

dU iW A 


: aero 1+ (y- If. [Pr - (2 





As for the stack region, integration is based on Eqs. (2.15), (2.17), and (2.19). In 
other words, the time averaged energy flux H, is only used in the stack region and is a 
constant throughout the stack region. 

The MATLAB numerical analysis method used in this work starts with specifying 
Re(p,), Im(p,), Re(U1), Im(U1), T_, Re(f), and Im(f) at a point in the duct. The complex 
pressure and volume velocity just outside the starting end of the ambient heat exchanger can 
be determined by integrating Eqs. (2.15) and (2.17) in the boundary layer approximation 


limit. Next, the complex pressure and volume velocity can be determined just inside the 


starting end of the ambient heat exchanger by invoking continuity of pressure and volume 
velocity at the end of the ambient heat exchanger. Then, Eqs. (2.17) and (2.18) are used to 
integrate over the ambient heat exchanger. Continuity of pressure and volume velocity is 
used to enter the stack. A value of H, is specified to get d7_/dx at the beginning of the 
stack. The solution in integrated to the heat exchanger. Again Eqs. (2.17) and (2.18) are 
integrated over the hot heat exchanger. Finally, the pressure, volume velocity, and the 


temperature at the exit of the hot heat exchanger are matched to those at the starting point. 


B. COMPLEX EIGENFREQUENCY 


The transition to onset in a prime mover is conveniently discussed in terms of the 
quality factor Q which can be defined as the ratio of 27 times the energy stored in the 
resonator to the energy dissipated per acoustic cycle. Letting W denote the dissipated 


power, 
=e (2.19) 


Both W and E,, are second order in the acoustic amplitude. A typical prime mover is 
comprised of five sections: the ambient duct, the ambient heat exchanger, the prime mover 
stack, the hot heat exchanger, and the hot duct. W is then the sum of the power dissipated 


in the five individual sections (Lin, 1989; Atchley, 1992 and 1994) 


W = W amb Si W ams hx + Wass a W ho hx + W tot ° (2.20) 


The subscripts amb, amb hx, pms, hot hx, and hot refer to the ambient duct, the ambient 


heat exchanger, the prime mover stack, the hot heat exchanger, and the hot duct, 
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respectively. The reader should refer to Atchley (Atchley et al., 1992) and Lin (Lin, 1989) 
for a full explanation. 

Instead of Q, we prefer to use the reciprocal of Q, because 1/Q converges to zero at 
onset rather than diverging to infinity. Also, it 1s this quantity that is proportional to the 
acoustic power output of the prime mover. It 1s worth mentioning that 1/Q is a function of 
the prime mover geometry, the thermophysical properties of the gas, p_, and the 
temperature difference. As the temperature difference along the stack increases from zero, 
the thermal losses in the stack decrease and eventually become negative, representing gain. 
In other words, below onset 1/Q is positive, at onset 1/Q is zero and above onset 1/O is 
negative. A complete evolution of Q through onset for a prime mover is discussed by 
Atchley (Atchley, 1994). 

In this work, the Q is determined by the complex eigenfrequency. Complex 
eigenvalues occur when a system has underdamped modes. The excitation of a resonator 
evolves in time as exp(i@,t) where the eigenfrequency is (Swift, 1988; Gamaletsos 1993; 
Arnott et al., 1994) 


@, = 2nf) + j Bo (22m 


Q 


where f, is the real resonance frequency. The quality factor Q thus can be found by 


Re(@, ) 


g= 2 Im(@, ) 


(2722) 


It should be noted that the complex eigenfrequency approach accounts for power 


generation in the stack and dissipation everywhere in the resonator. If this approach were 
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not used, an artificial acoustic source would have to be introduced into the system to deliver 


power to match boundary conditions (Choe, 1997). 


C. THE ANNULAR GEOMETRY FOR A PRIME MOVER 


As discussed in the previous chapter, the major difference between conventional prime 
mover geometry and that of an annular prime mover is that there are no typical dominating 
boundary conditions in the latter. Instead, an annular prime mover satisfies the periodic 


conditions 


P(r, 8) = p, (, +27), (2.23) 


and 


Ui(r, 0) = Ui(r, +22). (004) 


At first glance, this difference does not seem to cause much difficulty in the design and 
analysis of an annular prime mover, until one realizes the fact that almost all the previous 
work on conventional prime movers employ some dominating boundary condition as a 
convenient starting point of their analysis. Elimination of these boundary condition makes 
research on an annular prime mover considerably different from that of a conventional 
prime mover. 

It can be shown that an annular resonator can be modeled as a straight duct having an 
effect length L,, that is related to the inner and outer radii of the annulus and the particular 
mode under consideration (Choe, 1997). Therefore, we will treat the annular prime mover 
as a Straight duct having periodic boundary conditions. L,,.can be obtained by first writing 


the general solution of the wave equation in cylindrical coordinates which includes a linear 


Zh 


combination of the cylindrical Bessel and Neumann functions. Note that the Neumann 
function must be retained because the origin is excluded in the annular resonator. The 
effective circumference of the resonator can be then obtained by applying rigid boundary 
conditions at the inner and outer walls of the annulus. The annulus used in the experiments 
presented later has inner and outer radii of 10.0 and 15.3 cm, respectively. The height of 
the duct is 5.0cm. The cutoff frequency for cross modes is approximately 6.6 kHz. We 
are interested in the fundamental longitudinal mode (A = L,,) corresponding to frequencies 


in the 400 ~ 450 Hz range, well below cutoff. 


D. COMPUTER SIMULATION OF THE ANNULAR PRIME MOVER 


This subsection discusses the numerical analysis techniques for the annular prime 
mover. For the purpose of numerical integration in a thermoacoustic engine, Eqs. (2.7), 
(2.15), and (2.16) can be rewritten as a set of five, first order, ordinary differential 


equations in terms of five independent acoustic variables: Re(p,), Im(p,), Re(U1), Im(U1) 
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Our intention is to solve this set of coupled ordinary differential equations for Re(p;), 


Im(p,), Re U1), Im(U1) and T._,, subject to the periodic boundary conditions at the ends of 


the prime mover. This is the same basic method employed by DeltaE. 


1. Problem Statement 


When ordinary differential equations are required to satisfy boundary conditions at 
more than one value of the independent variable, the resulting problem is called a two point 
boundary value problem. Determining the deflection of a bar rigidly fixed at both ends is a 
typical example where the conditions specified are the deflections of the elastic curve at the 
supports. Heat-flow problems often fall into this class when temperature or temperature 
gradients are given at two points. The problem encountered here is of this kind, too. 

A major distinction between initial value problems and two point boundary value 
problems, as is pointed out by Press (Press et al., 1992), is that in the former case one is 
able to start an acceptable solution at its beginning (initial values) and just march it along by 
numerical integration to the end (final values). However in the latter case, the boundary 


conditions at the beginning point do not determine an unique solution to start with. For this 


og 


reason, two point boundary value problems require considerably more effort to solve than 
do initial value problems. Keller (Keller, 1992) has discussed the existence and 
uniqueness theory for two point boundary value problems. 

Two different numerical approaches to the annular prime mover are applied in this 


dissertation: DeltaE and a MATLAB program. 


2. DeltaE 


DeltaE is a computer program, developed by Bill Ward and Greg Swift at Los Almos 
National Laboratory, for modeling and designing thermoacoustic and other one- 
dimensional acoustic apparatus. Basically, DeltaE solves the one-dimensional wave 
equation based on the usual low-amplitude acoustic approximation. It numerically 
integrates the wave equation in a gas or fluid, in a geometry provided by the user as a 
sequence of segments, such as ducts, transducers, compliances, heat exchangers and 
stacks. It uses continuity of oscillating pressure, oscillating volume velocity, and mean 
temperature to pass from the end of one segment to the beginning of the next. It uses the 
appropriate wave equation and temperature equation for each segment. The iteration is 
controlled by global parameters such as frequency and mean pressure, and by local 
parameters such as the geometry of a segment and enthalpy flow. By defining a series of 
inputs, global parameters (the guess vector), the desired outputs, and boundary conditions 
(the target vector), the problem is solved iteratively. Guesses are updated until targets are 
achieved. The number of elements in the target vector must equal the number of elements 
in the guess vector, otherwise the system is over- or under-determined. When the 


computed and specified targets agree to within a specific tolerance, the solution is said to 
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converge. All calculations in DeltaE are performed in double precision and the user may set 
the error tolerance for matching targets. 

In general, because a pass of DeltaE’s integration does not reach the end with targeted 
values of all variables, a shooting method is used to adjust chosen initial variables in order 
to hit the target values. The user provides guesses for chosen initial variables. A 
successful convergence of DeltaE depends heavily upon a good choice of guess and target 
members. A good guide to the choice of guess and target variables is discussed by Ward 
and Swift (Ward and Swift, 1996). 

DeltaE is very versatile due to the fact that the elements of the guess vector are not 
limited to the conventional choices consisting of real and imaginary parts of p, and Ui. 
Any variables that have an effect on the target vector variables can be used as a guess. This 
feature enables DeltaE to compute resonance frequency, temperature, and geometrical 
dimensions in order to satisfy specified boundary conditions. DeltaE is very efficient at 
converging to solutions for some complicated acoustic systems, however, it knows nothing 
about acoustics or physics. For example, it does not recognize that negative frequencies, 
negative pressures or negative lengths are physically meaningless. It simply does the 
integration. Therefore, the reasonableness of the solutions produced will almost always 
depend on the quality of the initial guess vector and choices of the guess-target vector. 

To apply DeltaE to the annular prime mover the boundary conditions are incorporated 
into DeltaE’s target vector. The unknown conditions at the BEGIN segment, which DeltaE 
is supposed to find, are in DeltaE’s guess vector. To implement the periodic boundary 
conditions into DeltaE, a SOFTEND segment is inserted right after the initial BEGIN 
segment in the input file of the DeltaE program. Another SOFTEND segment (or they can 
both be HARDEND segments) is used at the end of the model. Four DIFFTARGETs are 


specified which match the amplitude and phase of p, and U; at these two SOFTEND 
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segments (Ward, 1997). The SOFTEND segments serve only as a tool to calculate the 
outputs and will be ignored in the target vector. SOFTENDs do not force the impedance to 
be zero. A given temperature across the prime mover stack is achieved by selecting the 
temperature at the hot heat exchanger as a target and the amount of heat input to the ambient 
heat exchanger as a guess vector. We now have five targets and only one guess; four more 
guesses are still required. The other four guesses have been selected tobe the amplitude 
and phase of p, and U, at the beginning segment. These guesses and targets (as are shown 


in Appendix B) are summarized in Table 2.1. 


1. Ip, (begin)! - Ip,(end)I < € 


5. Heat input at ambient heat exchanger 5. Hot heat exchanger temperature 











Table 2.1 Summary of the guess and target vectors in the DeltaE input file for the annular 


prime mover. € represents the tolerance for convergence. 


The important parameters for this work are the Q, the resonance frequency and the 
mode shape of the prime mover. To attain these desired quantities, more deliberation has to 
be made in determining the input file. Two methods can be utilized to find the Q with 
DeltaE (Ward, 1997). In the first method, the Q can be obtained indirectly by adding an 
artificial driver into the system, sweeping its frequency, and computing the frequency 


response of the system. The resonance frequency and half power points can be extracted 
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from this analysis. An alternate method is to compute the ratio of the stored energy to the 
dissipated energy. The stored energy is computed by adding an artificial driver, 
determining the pressure and velocity throughout the system and integrating the energy 
density. The dissipated energy is computed from the work done by the driver to excite the 
system. Since DeltaE only computes the pressure and velocity at the exit of a segment, the 
system must be broken into many segments to get a reasonable approximation to the stored 
energy. Thus the frequency response method was chosen to find the Q and resonance 
frequency. A side branch electroacoustic transducer was used as a driver, which has 
frequency independent parameters. A MATHEMATICA program (refer to Appendix A) 
was also developed to compute a least squares fit to the frequency response to find Q and 
resonance frequency. After finding resonance frequency, the system is then driven at 
resonance. To attain a detailed modal shape, the prime mover has to be divided into many 
small segments. The modal shape is a plot of the pressure at the end of each segments of 
the converged solution as a function of segment location. 

DeltaE is a very useful tool to predict how a given thermoacoustic apparatus will 
perform, or for helping the user design an apparatus to achieve a desired performance. 
However, it does have several disadvantages when applied to the annular problem. As was 
mentioned in Chapter I, in the presence of a stack, an annular prime mover exhibits 
frequency splitting. The degree to which a particular mode is excited depends on the 
location of the driver relative to the stack. Because of dissipation, Q for each mode is 
finite, which means that the two modes may overlap in the frequency domain. In terms of 
a driven acoustic system, this manifests itself as a multimode excitation. When only one 
mode is excited, it is straightforward to find Q and resonance frequency by a least squares 


fit to the frequency response curve. However, in the case when the driver excites two 


Zt 


mode simultaneously, some other technique, for example pole-zero analysis, must be used 
to determine Q and resonance frequency for each individual mode. 

Another disadvantage comes from trying to map out the mode shape with DeltaE. In 
DeltaE, values of acoustic pressure are only provided at the end of each segment. 
Therefore, to obtain a detailed mode shape, the prime mover must be divided into many 
small segments. This inherent property makes it is laborious and impractical to use DeltaE 
to find the mode shape. 

The other disadvantage is that currently there is no straightforward way to apply an 
externally-imposed temperature gradient on a resonator duct (Ward, 1997). This limitation 
makes it very difficult to model an annular prime mover and incorporate measured 
temperatures along the duct into the numerical model. 

Because of the disadvantages of applying the DeltaE program to our problem a 
decision was made to develop our own program. This program, written in MATLAB, was 
validated by comparing the results for a simplified annular prime mover with DeltaE 
results. 


3. MATLAB Program 


a. Approach 


One important aspect of this dissertation is to develop a numerical analysis 
program that is more tailored to the annular resonator problem than DeltaE. The program 
implements the shooting method because it appears to be a straightforward way to solve 
boundary value problems. This method also provides a systematic approach to taking a set 
of ranging shots that allows us to improve our aim systematically (Keller, 1992; Press et 


al., 1992). The program uses functions that already exist in MATLAB where possible. 
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A complete explanation of the procedure of the shooting method can be found in 
Press (1992) and Gerald (1994). Only a summary is presented here. Most differential 
equations of order higher than first can be reduced to coupled sets of first-order differential 
equations. The two point boundary value problem is equivalent to obtaining the solution to 
a set of N coupled first-order differential equations, satisfying N periodic boundary 
conditions at the starting point x, and at the end point x,. 


We start by writing Eqs. (2.25) through (2.29) in a general form 


dy;(x) 


rae Ce ee ea orn (2.30) 


with the periodic boundary conditions at the starting point x, and final point x, expressed as 


WX, Se eee), |. = Coen ey ees Vy | i= 12... Nes) 


where, in our problem, N equals five. 

To solve this boundary value problem, an initial value problem is created by 
assuming five initial values. The MATLAB function ode45, which implements fourth and 
fifth order Runge-Kutta formulas is then used to integrate the differential equations from x, 
to x,. The local error term for the fourth-order Runge-Kutta method is O(h’) and the global 
error would be O(h*), where h is the step size of the integration (Celia and Gary, 1992; 
Gerald and Wheatley 1994). Now, at the final point x,, a discrepancy vector F is defined 
with the dimension of N, whose components measure how far the first solution attempt is 
from satisfying the boundary conditions at x,. The problem now becomes one of finding 


the vector F such that 


Yi eee) Hoe nd lee, (2.32) 


ae 


To find the roots of this equation the Newton-Raphson method is employed. This method 
provides a very efficient means of converging to a solution, if a sufficiently good initial 
guess is given (Press, 1992). 

Let y denote the entire vector of values y,. In the neighborhood of y, each of the 


functions F, can be expanded in a Taylor series about y, as (Press, 1992) 


\. OF. 
Fy) = Flo +) = Flv) + LS* sy +0(G"). 2.33) 
j=l OY; 
The matrix of partial derivatives, known as the Jacobian matrix, J is defined as 
J,= =. (2.34) 
Since it is difficult to compute the matrix J analytically in this problem, finite differences 
are used to compute this matrix. In the MATLAB program, this is accomplished by 


perturbing each y, individually and finding the value of the vector of AF;/Ay,(x,). 


Equation (2.34) can be expressed in the matrix form 


F(y + dy) = Fly) + Jody +0(5”) . (2.35) 


To obtain a set of linear equations for the correction vector dy, terms of order dy’ 


and higher are neglected. Setting F(y+dy) = 0, gives 

Jody = - Fly). (2.36) 
Provided J is nonsingular, dy is given by 

dby = - J” F(y). (2.37) 
In practice, J'' is obtained by LU factorization. In MATLAB notation, dy is given by 
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dy = -J\F(y), (2.38) 


where \ is used to denote left division. 

The LU factorization implements Gaussian elimination. Because it involves the 
least amount of arithmetic operations, Gaussian elimination is generally considered to be 
one of the most efficient computation methods 1n the problem of solving a system of linear 
equations (Leon, 1994; Lindfield and Penny, 1995). It is worth pointing out that using 
J\F(y) instead of inv(J)*F(y) in MATLAB is two or three times faster and produces 
residuals on the order of machine accuracy relative to the magnitude of the data (The Math 


Works, 1994). Once dy is found, the solution vector is then updated as 


Ynew = Vo + OY. (2.39) 


The process is iterated until the solution converges, 1.e. when the value of the normalized 
norm of the discrepancy vector F is smaller than a desired threshold. 

In order to formulate a program to model the annular prime mover, the parameters 
to be updated at each pass of the integrations must be identified. These parameters are 
selected to be: the real and imaginary parts of the volume velocity U, and the complex 
frequency @, and the time-averaged energy flux along the stack H ,- Re(p,) and Im(p,) are 
kept constant at x = 0, which simply normalizes the amplitude of the pressure variation. 
In fact, for given values of Re(p,) and Im(p,), a unique solution can be obtained for a 
certain value of H,. The periodic boundary conditions are satisfied by matching the 
acoustic pressure p, and the acoustic volume velocity U, between the starting and final 
points. Also the temperature at the hot heat exchanger is matched to a certain value while 


keeping the temperature of the ambient heat exchanger a constant. 
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b. Description 


The complete MATLAB program is described in Appendix C. This section 
provides a brief explanation of how the program works. First, the coordinate system for 
the prime mover (as is shown in Fig. 2.2) is defined. The hot heat exchanger is at one end 
of the system. The end of the hot heat exchanger is atx = L,,,. The choice of the position 
of x = 0 is arbitrary, however by putting the stack and heat exchangers at the end, only 
four distinct regions are required. Had the stack been placed in the middle, there would 


have been five. 


Ambient Heat Prime Mover Hot Heat 
Exchanger Stack Exchanger 


Duct 











am, a ae 
1 Periodic = = = Periodic 
q eee 
‘ Boundary = = Boundary 
= SSS |S 
= SSE 

x=0 X=Legf 


Figure 2.2 Coordinate system used in the numerical analysis of the annular prime 
mover. 


Initial values of the real and imaginary parts of p,, a measured temperature profile 
T(x) along the duct and ambient heat exchanger, and a temperature for the hot heat 
exchanger, 7, are given. The real and imaginary parts of Ui, the real and imaginary parts 
of the resonance frequency @ and the time-averaged energy flux along the stack H, are 


guessed. It is assumed that the temperature of the duct and the ambient heat exchanger are 


BV 


held at certain values by external means. A set of reasonable initial guesses, which can be 
estimated from the output of DeltaE, is the key to successful convergence of the program. 
The desired temperature of the hot heat exchanger (7, in the program) is accomplished by 
adjusting the time-averaged energy flux along the stack H, . The targets to be matched are 
the real and imaginary parts of p, and U, at x = 0 and x = Ls, and the temperature of the 
hot heat exchanger. This matching is accomplished by adjusting’ the real and imaginary 
parts of U1, the real and imaginary parts the complex angular frequency @, and the time- 
averaged energy flux along the stack. It is also assumed that there are no temperature 
gradients in the ambient and hot heat exchanger. 

Once the required input parameters are determined. The program solves the 
system of equations for the five unknown variables. If the resulting new p,, U1, and T, 
differ significantly from the initial guesses it is necessary to update the guesses and do the 
integration again. During the integration, all the temperature dependent parameters and gas 
parameters are calculated according to the temperature distribution. This whole process is 


repeated to convergence. 
c. Validation 


This program was validated by comparing its solution to that using DeltaE for a 
simple case, in which it is assumed that the temperatures of the entire duct and the ambient 
heat exchanger are held at 293 K by some external means. The desired temperature of the 
hot heat exchanger is achieved by adjusting the quantity H ag 

The comparisons of the results from the two programs are shown ‘in Figs. 2.3 
through 2.8. Figure 2.3 shows resonance frequency vs. AT, with AT increasing from 


OK to 240 K. Figure 2.4 shows 1/Q vs. AT over the same temperature difference span 


315, 


as in Fig. 2.3. The overall agreement of resonance frequencies and 1/Q from the two 
methods is good. In fact, the normalized mean square errors of the resonance frequency 
between two methods are 1.48x10°% and 6.10x10°% for the low and high frequency 
mode, respectively. The normalized mean square error of 1/Q between two methods are 
0.684% (low frequency mode) and 0.854%(high frequency mode). The value of 1/O 
calculated from DeltaE is determined from a pole-zero analysis of the frequency response. 
The value of Q so determined is sensitive to the bandwidth and weighting used in the fit. 
The error bar attached to the AT = 200 K data point indicates the uncertainty and accounts 
for most of the disagreement at higher values of AT. One reason for the dependence of QO 
on bandwidth is the multimode excitation. A driver was incorporated in the DeltaE 
program to find the frequency response, from which we obtain the required modal 
information. The driver locations are selected such that at AT = 0K _ there is only one 
mode being excited. However, as AT is increased, the mode shapes can be altered and the 
previous driver position may excite both modes simultaneously. This can be seen in Figs. 
2.5 and 2.6 which show the frequency response of the low mode and high mode of the 
prime mover from DeltaE for three values of AT. For higher the AT, the effect of high 
mode excitation is clearly visible as a small bump in the frequency response of the low 
mode. However, no such bump is observed in the frequency response of the high mode. 
Figures 2.7 and 2.8 show the mode shapes for the low frequency mode with 
AT=OK and AT=100K. _ Figures 2.9 and 2.10 show the mode shapes for high 
frequency mode with AT=0K and AT=100K. In these figures, the solid lines 
represent the results of the MATLAB program and the open circles represent the results of 
the DeltaE program. The mode shapes from DeltaE have been normalized to the highest 
value of the pressure amplitude. The amplitude of the mode shape from the MATLAB 


program is determined from a least squares fit to the DeltaE results. 
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It is seen that one of the advantages of the MATLAB program is the smooth 
mode shape curve. The agreement of the mode shapes are good in general. However the 
agreement of the low frequency mode is not as good as for the high frequency mode, 
particularly at high temperature difference. The reasons for this discrepancy at higher 
values of AT are not fully understood. It should be pointed out that the two programs treat 
the heat exchangers differently. The MATLAB program assumes that the mean gas 
temperature at the heat exchanger equals the temperature of the heat exchanger. However, 
DeltaE imposes a thermal impedance so that the gas and heat exchanger temperatures differ. 

Figure 2.8 displays one interesting feature. Both program predict a low standing 
wave ratio. This feature is not present in the high mode. 

Although there remain small disagreements between DeltaE and the MATLAB 
program, the major features are in good agreement. In the chapters to follow, the 
MATLAB program will be used to analyze the annular prime mover under more realistic 
conditions. The major difference is that the measured temperature profiles are incorporated 


in the program. 
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Calculated resonance frequencies VS DeltaT for the annular prime mover 
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Figure 2.3 Comparison of the calculated resonance frequencies of the annular prime 
mover as determined from the MATLAB program and DeltaE. Symbols are explained 


in the legend. Agreement 1s excellent at all values of AT. The normalized mean square 


error is 1.48x10°% for the low mode and 6.10x10°% for the high mode. 
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Calculated 1/Q VS DeltaT for the annular prime mover 
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Figure 2.4 Comparison of 1/Q of the prime mover vs. AT as determined from the 


MATLAB and DeltaE. An error bar is indicated for the low mode at AT = 200 K. 
Because of modal overlap, determination of 1/Q in DeltaE is not without uncertainty. 
Agreement is better at low AT. The normalized mean square error is 0.684% for the 


low mode and 0.854% for the high mode. 
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Frequency response of the low mode for the annular prime mover 
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Figure 2.5 Frequency response of the prime mover from DeltaE (low mode) at 


AT=0OK, 100K and 200K. Note the small bump at about 436 Hz on 


AT = 100 K and 200K curves which is a result of mode overlap. 
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Frequency response of the high mode for the annular prime mover 
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Figure 2.6 Frequency response of the prime mover from DeltaE (high mode) at 


AT =0 K, 100 K and 200 K. Notice that there is no visible bump from overlap. 
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Low mode with DeltaT = 0 Kelvin ; frequency = 419.2+6.038i 
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Figure 2.7 Comparison of mode shapes of the prime mover at AT = 0 K_ for the low 


frequency mode as calculated with the MATLAB program and DeltaE. The beginning 
of the stack region is indicated by the dashed vertical line. 
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Low mode with DeltaT = 100 Kelvin ; frequency = 422.3+6.192i 
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Figure 2.8 Comparison of mode shapes of the prime mover at AT = 100 K for the 


low frequency mode as calculated with the MATLAB program and DeltaE. Although 


both programs predict the same general mode shape, the agreement is not as good at 


AT=0OK. 
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High mode with dT = 0 Kelvin ; freq = 437.1+2.013i 
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Figure 2.9 Comparison of mode shapes of the prime mover at AT = 0 K for the high 


frequency mode as calculated with the MATLAB program and DeltaE. 
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High mode with dT = 100 Kelvin ; freq = 4374+2.137i 
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Figure 2.10 Comparison of mode shapes of the prime mover at AT = 100 K for the 


high frequency mode as calculated with the MATLAB program and DeltaE. The 


agreement is much better than for the low mode at AT = 100 K. 
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EK. END EFFECTS 


The usual boundary conditions applied to a change in cross section are continuity of 
pressure and volume velocity. However, at the discontinuity an additional acoustic 
impedance is introduced owing to the abrupt change of cross sectional area (Miles, 1946 
and 1947; Karal, 1953; Pierce, 1991). As a matter of fact, we have shown (Choe, 1997) 
that it is critical to take into account end corrections at the constriction to get good 
agreement between the measured and predicted resonance frequency and mode shapes for a 
constricted annular resonator. 

Sudden changes in duct cross section produce sudden changes in acoustic pressure 
amplitude. Below cutoff, the effect 1s equivalent to a lumped-parameter impedance. For 
example, a constriction in a duct gives nse to a concentration of flow through the 
constriction that consequently increases the kinetic energy density. As is discussed by 
Morse and Ingard (1968), if the net increase is concentrated in a region small compared to a 
wavelength, the total increase can be characterized by a lumped acoustic inductance. If the 
constriction also produces an increase in viscous energy loss, this addition can be modeled 
as a lumped acoustic resistance. 

This equivalent lumped impedance can be obtained from various methods. Only two 
methods, the conformal transformation method and the higher order mode method, are 
discussed here. 

After computing the analogous acoustical impedance Z, the boundary condition for the 


pressure becomes. 


p.() = p,() - Z*U,(), (2.40) 
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where p, is the pressure in the unconstricted region and p, is the pressure in the 
constriction. U, is the volume velocity in the unconstricted region and the constriction 


junction is at x =/. The volume velocity boundary condition is unchanged. 


1. Conformal Transformation Method 


The first method used to compute the lumped impedance introduced by a constriction 
is discussed by Morse and Ingard (1968). To obtain the analogous impedance of a 
constriction, they compute the extra kinetic energy and viscous-energy loss produced by 
steady flow through the constriction. To find the steady-state flow, the interior of the half- 
duct is first transformed to the upper half of a complex w plane by the Schwartz-Christoffel 
transformation. Next, depending on the geometry of a specific problem, a further 
transformation is required to go from the w plane to the velocity-potential plane in which 
the lumped impedance can be readily obtained. 

Morse and Ingard applied this technique to the case of a rectangular duct of depth d 
and width b for x <0 and a for x >O (as shown in Fig. 2.11). The lumped-circuit 
inductance and resistance corresponding to the sudden increase in the cross section of a 


rectangular duct are 














2 a 
— b 
L, = hl (a-b) | a+b ui 1p ete) (2.41) 
td| 2ab a-—b 4ab 
and 
2 ee 
pe =O AE 5 Gs Pa aea (2.42) 
2ad_ b Ttab a-b | 


Notice that both L, and R, vanish when a = b (1.e., when there is no change in width). 
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Figure 2.11 Section of a rectangular duct with sudden change of width in the y 
direction. 


2. Higher Order Mode Method 


This method involves including higher order modes near the constriction ends 
(Muehleisen, 1996). In this work, Muehleisen examines the reflection, transmission, 
radiation, and coupling of higher order modes at a discontinuity in finite length rigid wall 
rectangular ducts. Using a method of generalized scattering coefficients, an analytic 
expression for the reflection and transmission of higher modes at size discontinuities is 
developed. 

From the scattering matrix formulation an expression for the lumped acoustic 
impedance can be found (Muehleisen, 1997). This method does not include losses, 


however, and so the resistance obtained by conformal mapping is added to the impedance. 
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WT. EXPERIMENTAL APPARATUS AND PROCEDURE 


A. INTRODUCTION 


Descriptions of the experimental apparatus and data acquisition procedure are 
presented in this chapter. The chapter begins with a discussion of the annular resonator. 
This resonator was used in preliminary work to measure the modal properties of a 
constricted annular resonator in the absence of a stack (Choe, 1997). This discussion is 
followed by a detailed description of the annular prime mover. Next, the constricted 
annular prime mover is discussed, followed by a description of how the prime mover is 
instrumented with microphones and thermocouples. The calibration of the microphones is 
discussed in this subsection as well. The chapter concludes with a description of the 


procedure followed during data acquisition. 


B. EXPERIMENTAL APPARATUS 
1. Annular Resonator System 


The annular resonator, shown in Fig. 3.1 (without the top) and Fig. 3.2 (with the top) 
has been used in a preliminary work to investigate the modal properties of a constricted 
annular resonator (Choe, 1997). The inner and outer radius of the resonator is 
10.00 + 0.01 cm and 15.30 + 0.01 cm, respectively. The outer side wall is a 1 cm thick 
aluminum ring with 15.30 + 0.01 cm inner radius. A solid aluminum 10cm _ radius 
cylinder that is concentric with the outer ring forms the inner side wall. The height of the 


resonator is 5.00 + 0.01 cm. The dimensions of the resonator were chosen to give a 
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longitudinal resonance frequency in the 300 to 500 Hz range, which is well below the 
cutoff frequency for the first cross mode, as discussed in Chapter II. 

Since the prime mover will be operated at temperatures over 200 °C, the plexi-glass 
top described in Choe has been replaced by a 0.64 + 0.01 cm thick aluminum top with a 
radius of 16.01 + 0.01 cm. The top is bolted with a 10.7 cm long threaded rod at the 
center and three equally-spaced 7.1 cm long clamp bars on the edges. Ten microphones 
are installed in the top to measure the spatial pressure distribution inside the driven 
resonator below onset and observe the behavior of the sound generated above onset. To 
monitor the temperatures in the resonator, 5 type-E thermocouples are installed through the 
top and another 5 type-E thermocouples are installed in the stack/heat exchanger assembly. 
A detailed description of the microphone and thermocouple installation will be presented 
later in this chapter. 

It is very difficult to move the stack and heat exchanger assembly after it has been put 
into the resonator because of the complex wiring and the fragility of the glass components. 
Therefore, three holes (3.50 mm diameter) were drilled in the bottom of the resonator to 
allow the driver location to be moved relative to the stack. Depending on which mode is to 
be excited, the driver was then connected to the desired hole viaa0.4cm O.D., 3.5 cm 
long plastic tube. The two unused holes were plugged with thumbtacks that had some 
Dow Corning high vacuum grease applied to the tips. 

To ensure tight sealing for a high Q, a thin layer of Dow Corning high vacuum grease 
was applied to the surface of the aluminum ring and the center piece before the top was put 
on and bolted down. In addition to the bolt in the center and the three clamps mentioned 
earlier, 6 extra equally spaced clamp bars are applied around the edge of the top to ensure a 


tight seal. 
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The instrumentation is shown in Fig. 3.3. The acoustic signal is generated by a 40W, 
16 €2 external compression driver. The output of a Hewlett Packard 3562A dynamic signal 
analyzer was sent through a Techron 5507 power amplifier to the driver. The signal 
analyzer was used to determine the frequency response of the resonator. Typically, the 
resonator was driven with a random noise source of 1 Vrms source level. The pole-zero 
analysis function built into the signal analyzer was then used to determine the Q and the 
resonance frequency. To measure the mode shape, the source of the analyzer was changed 
to a 1Vrms fixed sinusoidal signal at the resonance frequency found in the previous 


measurement. 


2. Stack Assembly 


The prime mover stack assembly consists of the stack plates and the stack holder. The 
first challenge was to select a suitable material for stack. The stack needed to withstand 
high temperatures and the differential expansion imposed upon it from the high gradients. 
Glass has a low thermal conductivity in comparison to most metals (Johanson Companies, 
1996) and is able to withstand higher temperatures than most plastics and other polymers. 
An AF-45 thin borosilicate glass slide (made by Precision Glass & Optics) was chosen to 
be the material for stack. Some properties of this AF-45 glass are provided in Appendix D. 
AF-45 is a modified borosilicate glass with a high content of BaO and AlI,O, and features a 
low thermal conductivity and a low coefficient of thermal expansion of 45.0 x 107K"! 
(20~300° C). The glass was accurately cut to a size of: 43.65+/-0.05 mm wide and 
23.90+/-0.10 mm long and is 0.145+/-0.02 mm thick. 

The stack holder (detailed dimensions of which are given in Appendix E) consists of 


two horizontal bars (5.28 cm long) connected to two vertical bars (4.85 cm long) as shown 
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in Fig. 3.4. Tongue-and-groove joints are machined in the ends of the bars to ensure that 
they form a sturdy frame. Sixty-one 0.508 mm deep, 0.203 mm wide, equally-spaced slits 
were cut in the inside of each of the two horizontal bars. The center-to-center spacing 
between the slits is 0.762 mm. The gaps between the cover glass plates are 0.610 mm. 
After the stack holder was assembled, 61 glass slides were then carefully inserted into the 
slits. The overall porosity for the prime mover stack assembly, including the frame, is 


approximately 0.61. 


3. Ambient Heat Exchanger Assembly 


The ambient heat exchanger is used to remove heat from the stack and to maintain one 
end of the stack close to ambient temperature. The ambient heat exchanger assembly is 
constructed as shown in Fig. 3.5 and the detailed dimensions are given in Appendix E. 
The heat exchanger 1s designed to ensure that the temperature difference between the center 
and the edge of the heat exchanger is less than one degree C fora 35 W load. The ambient 
heat exchanger holder was made in the same way as the stack holder, except that it uses 
copper instead of aluminum and has only 10 slits. Also, a central bar was incorporated to 
supplement the heat transfer between the heat exchanger and the resonator body. The slits 
are 0.635 mm deep, 0.711 mm wide, and equally spaced at 4.064 mm (center to center). 
After the holder was made, twenty Alloy-100, 2.30 mm long, 5.10 mm wide, and 0.68 
mm thick copper fins were cut and soft-soldered into the slits (10 on either side of the 
central bar). The overall porosity for the ambient heat exchanger is 0.62 (including the 
frame). The next step was to soft-solder a 5.256 cm x 4.977 cm rectangular copper woven 


wire cloth (20 x 20 mesh) on one side of the heat exchanger assembly. The function of 
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this wire cloth is to increase the thermal contact between the prime mover stack and the 
ambient heat exchanger. 

The fin spacing is much larger than would be used in a practical prime mover. The 
acoustical properties of an constricted annular resonator are very sensitive to both the 
porosity and flow impedance of the constriction. It was decided to make the flow 
impedance of the heat exchangers small compared to that of the stack. This choice 


necessitated using large fin spacing, which reduced the effectiveness of the heat exchanger. 


4. Hot Heat Exchanger Assembly 


The function of the hot heat exchanger is to supply heat to the stack. It consists of 
nickel-chromium wire woven on a non-conducting frame. Since nickel-chromium wire is 
used to provide heat to the system, the hot heat exchanger has to provide electrical 
insulation to the aluminum resonator. Because of the need for electrical insulation and high 
temperature operation, virgin electrical grade PTFE (Teflon) is chosen for the hot heat 
exchanger frame. The hot heat exchanger assembly is shown in Fig. 3.6. Drawings of 
fabricated parts are provided in Appendix E. The Teflon frame has outer dimensions of 
4.978 cm wide, 5.283 cm high, and 0.813 cm thick. The inside of this frame was milled 
out to a size of 4.470 cm wide and 4.663 cm high. Seventeen holes are drilled onto the 
two vertical sides with an equal spacing of 0.521 cm (8 on one side and 9 on the other 
side). Into each of the holes is inserted a 0-80 stainless steel screw. A 92.1 cm, # 28 gage 
nickel-chromium wire with a total resistance of 12 42 was wound around the 17, 0-80 
screws. To avoid direct contact between the nickel-chromium wire and the stack, stainless 
steel flat washers are used on the screws for spacers. The ends of the nickel-chromium 


wire are wrapped with heat shrinkable flexible PVC tubing to insulate the wire from the 


2 


resonator. The wire is then fed through Teflon tubing which is placed and epoxied in the 


bottom wall of the resonator. 
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Figure 3.1 The annular resonator (w 
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Figure 3.2 The annular resonator (w 
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Figure 3.3 Block Diagram of Instrumentation Setup. 
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Figure 3.5 The ambient heat exchanger assembly. 
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Figure 3.6 The hot heat exchanger assembly. 
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5. Stack/Heat Exchanger Assembly 


The stack assembly, the ambient heat exchanger assembly, and the hot heat exchanger 
assembly were connected with four 5.08 cm long, 00-90 stainless steel threaded rods on 
the two horizontal sides (as shown in Fig. 3.7). Teflon spacers were used to ensure proper 
spacing between the stack and the heat exchangers. The spacing between the glass plate 
and the nickel-chromium wire is approximate | mm. Before the stack/heat exchanger 
assembly was placed into the resonator, a piece of thin Teflon tape was wrapped on the 
side of the hot heat exchanger to avoid electrical short circuit. A thin layer of silicone heat 
sink compound was applied to the side of the ambient heat exchanger frame to supplement 
the heat transfer between the ambient heat exchanger and the resonator wall and to seal any 


gaps between the frame and the resonator. 
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Figure 3.7 The stack/heat exchanger assembly. 
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6. Constriction 


Placing another constriction in the resonator provides the possibility for an annular 
prime mover to reach onset. Such constrictions were previously used in related work by 
Choe (Choe, 1997). The constrictions are made of PVC (Ployvinylchloride) because of its 
rigidity and ease of fabrication. The area ratio is defined by the ratio of open area in the 
constriction (which 1s the shaded block in Fig. 3.8) to the cross-sectional area of the 
resonator. In Fig. 3.8, S, stands for the open area and $, stands for the open cross- 
sectional area of the constriction. The table in Figure 3.8 gives the geometrical parameters 
of the constrictions. Constrictions were used having area ratios of 0.7, 0.3 and 0.1 with a 
total length of 45°. The 45° length was achieved by stacking two 11.25° constrictions and 
one 22.5° constriction together. Vacuum grease was applied to the joining sides of the 


individual constrictions to ensure a tight fit. 
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Figure 3.8(a) Nominal geometry of the constrictions. 
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Figure 3.8(b) Measured geometry of the constrictions. 
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7. Microphones 


In this experiment, 8 Panasonic electret WM-60AT microphones and two ENDEVCO 
piezoresistive pressure transducers (Model #8510B-5, Serial # G63T, sensitivity 53.85 
mV/psi, and Model #8530C-15, Serial # AGLC2, sensitivity 10.89 mV/psi) are used to 
measure the pressure distribution inside the resonator. The two ENDEVCO piezoresistive 
pressure transducers are intended to be used to observe the behavior of the sound generated 


above onset, in case the Pansonic microphones saturate. The Panasonic microphones are 6 
mm in height and 5 mm in diameter, with a sensitivity of -44 dB + 3 dB re 1 V/Pa. A 
Hewlett Packard 6237 Triple Output Power Supply provided 9V DC to each Panasonic 
microphone via a custom made microphone selector. 


As shown in Fig. 3.9, the 8 Panasonic microphones were spaced by 45 degrees in 


pre-drilled holes (0.605 cm in diameter, 0.50 cm deep). They are glued into the holes with 
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silicone adhesive. At the bottoms of these holes, 0.82 mm diameter holes were drilled to 
expose the microphones to the resonator cavity. The two ENDEVCO piezoresistive 
pressure transducers were flush mounted in the cap through spacers with 0.404 cm in 


diameter threaded holes in the center. 


Qdegrees 45 degrees 







PANASONIC 
Microphones 


mic #7 





90 degrees 


mic #6 


ENDEVCO Pressure 


Transducers 
180 degrees 


Figure 3.9 Installation of 10 microphones on the resonator top. 


Prior to the experiment, the 8 Panasonic microphones were calibrated relative to 
Microphone # 1. The calibration apparatus is shown in Fig. 3.10. A calibrating chamber 
was made from a 19.6 mm long, 21.55 mm ID, and 27.22 mm OD circular PVC pipe. A 
Panasonic micro-speaker (Model #EAS-2P106C) was epoxied into one end of the pipe. 
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The resonator top was removed and the calibration chamber was positioned over each 
microphone. A thin layer of Dow Corning high vacuum grease 1s applied on the edge of 
the chamber to ensure a tight sealing. The frequency responses were then measured for 
each microphone over a 300 ~ 500 Hz frequency range using a Hewlett Packard 35665A 
Dynamic Signal Analyzer. For the calibration, the signal analyzer was set to swept sine 
mode with a integration time of 20 cycles, a settling time of 20 cycles, and a resolution of 


200 pts/sweep. 


PANASONIC 
Micro-Speaker 











Calibrating Chamber | HP 35665A DYNAMIC 
PANASONIC 
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Figure 3.10 Calibration of the Panasonic microphones using a small chamber. 


Microphone #1 was selected as the reference for the comparison calibration. The 
frequency response of microphone #1 when covered by the calibration chamber is shown 
in Fig. 3.11. Each microphone’s frequency response was divided by that of microphone 


#1 to create a correction factor. A typical curve for the correction factor is shown in Fig. 
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3.12. After finishing the calibrations, one more measurement was made on Microphone #1 


to ensure the reproducibility of the data. With this particular calibration chamber, a 
difference of less than 1% was achieved. The correction factor for the 8 microphones is 


usually in the range of 0.75 to 1.75. 

When measuring the mode shape, the output from each microphone was divided by 
the appropriate correction factor to obtain the corrected amplitude at a particular location and 
driving frequency. A custom made microphone switch board selects each microphone to 
measure the pressure amplitude at different locations. The microphone output amplitude is 


measured using a Standford Research 530 Lock-in amplifier. 


 X:400 Hz NO SVB 
Freq Resp 2:1 





300 Hz 350 Hz 


Figure 3.11 Frequency response of microphone #1 when covered by the small 
calibration chamber. 
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© X:400 Hz Y¥:939.5697 m 
FRES/D1 2:1 





Figure 3.12 Plot of correction factor of microphone #2 referenced to microphone #1. 
The correction factor at 400 Hz is 0.940. 


8. Thermocouple Instrumentation 


Ten Omega type E (Model # S5TC-GC-E-30-72) chromel-constantan 0.254 mm 
diameter glass braid insulated thermocouples are used to monitor the temperatures of the 
Stack, the heat exchanger, and inside the resonator. Type E thermocouples were selected 


over other types because they provide the greatest sensitivity (in V/°C) and because of their 


useful temperature range. 

As shown in Fig. 3.13, five thermocouples are used to monitor the temperature of the 
stack and the ambient heat exchanger. Two thermocouples each are placed on the hot and 
ambient sides of the center glass slide, one in the center and the other at the edge. These 


thermocouples are fed through 1.4 mm diameter holes on the copper woven wire cloth of 
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the ambient heat exchanger. Another thermocouple (thermocouple # 6 in Fig. 3.13) was 
placed on one of the fins of the ambient heat exchanger. The thermocouples were attached 
to the glass plate using JB weld adhesive because its high temperature capability. The size 
of the glue spots is approximate 0.01 mm’. These glue spots are approximate 0.5 mm 
from the edge the glass plate. The wires of the five thermocouples were individually 
threaded through 1.4 mm diameter holes in the resonator bottom adjacent to the ambient 
heat exchanger. The other five thermocouples were used to sense temperature inside the 
resonator, as is shown in Fig. 3.14. They were fed through 1.4 mm diameter holes in the 
top and positioned approximately 3 cm inside the resonator. All the holes in the resonator 
were then filled with BIPAX TRA-BOND epoxy to ensure a tight seal. 

The temperature difference across the stack is adjusted by the voltage supplied to the 
hot heat exchanger which is controlled by a Power Design 3650R DC power supply. The 
temperatures inside the resonator and in the stack are measured by a Keithley 740 system 


scanning thermometer. 
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Figure 3.13 Thermocouple placement on the stack and the ambient heat exchanger. 
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Figure 3.14 Thermocouple placement in the resonator top. 


C. EXPERIMENTAL PROCEDURE 
1. Resonator Setup 


The stack/heat exchanger assembly is placed at a fixed location relative to the 
microphones and thermocouples. For the experiment, three driver positions are used: 45°, 
90°, and 180° from the center line of the stack/heat exchanger assembly. Untlted driver 
holes are plugged with thumbtacks with some vacuum grease applied to the tips. The 


resonator top is oriented so that 0° is set at the stack assembly center line. The O° setting is 
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performed by rotation of the top by hand. To determine the error in positioning the top, the 
top was positioned five times. It was found that the error by rotation was less than + 1.0°. 
So, with the mechanical error + 0.25° in drilling the microphone holes, the total error of 
microphone positioning is less than + 1.25°. With the stack/heat exchanger assembly 
positioned, the proper amount of vacuum grease is applied on the open surface of the 
resonator and the top is secured with the clamp bars. 

For the measurements with a constriction in the prime mover, a 45° long constriction 
is inserted at a location that is 90° from the center line of the stack/heat exchanger assembly 


to the center of constriction. 


2. Determining Resonator Characteristics 


After the resonator was set, the Hewlett Packard 3562A Dynamic Signal Analyzer is 
used to characterize the resonator. The swept sine mode of the analyzer provides a sine 
wave signal to the system at a given range of frequencies. The frequency response of the 
system is used to determine the resonance frequency and Q of several modes of the 
resonator. When only one mode is excited, the peak pressure amplitude, frequency, and 
quality factor are same as the those given by pole-zero analysis. When two modes are 
excited simultaneously, pole-zero analysis is used to extract the resonance frequency and 
quality factor. The signal analyzer gives poles in the form a+) b. The resonance is simply 
the value of b and the quality factor can be obtained by b/2a. This method has been 
explained by Choe (1997). These measurements are performed for (1) annular prime 
mover without a constriction, and (2) annular prime mover with three different 
constrictions (45° long, porosities of 0.7, 0.3 and 0.1) placed at 90° relative to the center 


of the stack/heat exchanger assembly. The driver in the first experiment was located at 90° 
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(to excite the low frequency mode) or 180° (to excite the high frequency mode) from the 
center of the stack assembly. In experiment (2), the driver was placed at 45° from the 


center of the stack assembly, between the stack and constriction. 


3. Data Collection 


Data collection 1s performed in sets of annular prime mover -without and with a 
constriction. Before data collection, the output of the Hewlett Packard 6237B Triple output 
power supply is set to DC 9V. The output voltage of the Power Design 3650R power 
supply is also adjusted to a certain value and the temperature at the hot end of the stack is 
allowed to reach steady state (temperature variation < 0.1°C/10 sec). The measurement is 
commenced by recording the temperatures measured by the Keithley 740 System Scanning 
thermometer. The resonance frequency and quality factor is measured by the Hewlett 
Packward 3562A Dynamic Signal Analyzer and the microphone switch board is set to the 
microphone that provides the maximum output. The procedure used to acquire the 
resonance frequency and quality factor was described in subsection 1. After finding the 
resonance frequency, a data run is commenced by changing the source of the Hewlett 
Packward 3562A Dynamic Signal Analyzer to a 1 V,, fixed sine with a driving frequency 
equal to the resonance frequency obtained by previous procedure. The acoustic amplitude 
iS Measured with each microphone. With the driver running, the switch board selects a 
particular microphone and the value of amplitude at each position is measured with the 
Stanford SR530 Lock-In amplifier. The Lock-In amplifier is set to magnitude-phase 
display mode with a time constant of 1 second. At this time, the temperatures are recorded 
again and same measurement are repeated again. Before the entire set of experiments is 


finished, the final temperatures are recorded. Then, the output voltage of the Power Design 
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3650R power supply is changed to adjust the hot end temperature and the same procedure 
is followed. All data recorded from the 8 microphones are scaled by the correction factor 
obtained from the microphone calibration. The corrected data are then entered in a 


computer program to plot the mode shape. This concludes data acquisition of one mode. 


WZ 


IV. RESULTS AND DISCUSSION 


This chapter shows the comparison between theoretical and measured results for two 
series Of experiments on the annular prime mover. First, measured temperatures in the 
prime mover are presented. Next, results are presented for resonance frequency, 1/Q vs. 
AT, and mode shape for the unconstricted annular prime mover. Results for the annular 
prime mover with a constriction located 90° from the center of the stack/heat exchanger 
assembly are then discussed. Finally, an analysis of the potential for onset in a two-stack 


annular prime mover is presented. 


A. ANNULAR PRIME MOVER 


The measured resonance frequencies and Q’s are determined as functions of AT from 
pole-zero analysis of the frequency response. The theoretical resonance frequencies and 
quality factors are determined from the complex frequency calculated using the measured 
temperature distribution in the prime mover as an input. Figures 4.1 and 4.2 show the 
measured temperatures as functions of the heater voltage at four locations on the center 
plate of the stack for the low mode and high mode of the prime mover. 

Figure 4.1 shows that the transverse (center-to-edge) temperature difference on the 
ambient side of the stack is at most 5.4 K. However, there is a larger transverse 
temperature profile on the hot side. At the highest heater voltage the hot-side, center-to- 
edge, temperature difference is 14.8 K. The transverse temperature profile leads to some 
uncertainty in the value used for AT. The effects of a transverse AT profile on the 


performance of a stack is difficult to access. Therefore, the transverse profiles will be used 


We 


to place error bars on the AT values. Figure 4.2 shows the center-to-edge temperature 
difference for the high mode. It is seen that the ambient-side temperatures are essentially 
the same as in Fig. 4.1. However the hot side temperatures show a larger center-to-edge 
difference than in Fig. 4.1. The difference is attributed to the fact that the resonator top had 
to be removed to reposition the driver to switch from exciting the low mode to the high 
mode. It is entirely possible, and observed, that while removing and replacing the top the 
stack/heat exchanger assembly was disturbed. This could have slightly altered the hot heat 
exchanger/stack spacing, resulting in the observed differences. 

In the results to follow, all the low mode data were taken, then the driver was 
repositioned and all the high mode data taken. Therefore, the temperatures displayed in 
Figs. 4.1 and 4.2 should be representative. Typical measured temperature distributions in 
the prime mover are shown in Figs. 4.3 (the low mode) and 4.4 (the high mode). The 
measured temperature profile in the duct (0 to 0.792 m) and the ambient heat exchanger 
temperature are used as inputs into the MATLAB program. The measured hot end stack 
temperature is a target. 

Figure 4.5 shows the comparison of the measured and computed resonance frequency 
for the low mode and high mode. The driver location is 90° from the center of the 
stack/heat exchanger assembly for the low frequency mode and 180° for the high frequency 
cle It is seen that there is good agreement between the measured and predicted 
resonance frequencies. The biggest differences between the measured and predicted 
resonance frequencies are approximately 0.86 % and 1.77 % for the low mode and high 
mode, respectively. 

Figure 4.6 shows 1/Q of the low and high mode for the annular prime mover verses 
the temperature differences across the stack. In general the agreement is good, although 


there is a tendency to under predict the 1/Q (or over predict Q), indicating the presence of 
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unaccounted losses. For the low mode, 1/Q decreases as the temperature difference is 
increased while for the high mode 1/Q increases as the temperature increases. It is evident 
from these results that the annular prime mover will not reach onset under readily attainable 
conditions. As a matter of fact, an extrapolation of the 1/Q curve for the low mode gives an 
onset temperature difference of approximately 1400 K. 

The larger discrepancies of the resonance frequency and-1/@ for the high mode at high 
temperature gradient can be caused by a misjudging of the stack temperature. In the 
original calculations, the temperatures of the hot (7,) and ambient (7,) sides of the prime 
mover stack are calculated by a weighted average of the, temperatures at the edge and center 


of either side of the stack 


Zit 1. 


T = ——e he 4,1 
h g ( ) 
and 
27 1 
T= ne - =, (4.2) 


where 7,. and 7,. are the center and edge temperature of the hot side of the stack, 
respectively. 7,. and 7,, are the center and edge temperature of the ambient side of the stack, 
respectively. If the simple average of 7, and T, are used in the computations as opposed to 
the weighted average, the resonance frequencies are almost unchanged (difference < 
0.03%). This negligible change is to be expected since the stack/heat exchanger assembly 
is very short compared to the effective length of the prime mover 
(Lesack assembly ! Leg = 0-038). As a results, the vertical error bars in Fig 4.5 are smaller 
than the symbols used to plot the data. However, 1/Q is more sensitive to uncertainties in 


T, and 7, than is the resonance frequency since all thermoviscous effects are most 


important in the stack region. The sensitivity to uncertainties in AT is represented by the 
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vertical error bars shown in Fig. 4.6. There is also an uncertainty in the calculation of 
temperature difference across the stack. As shown in Fig. 4.1 and 4.2, the temperatures at 
the center of the hot end and ambient end differ from those at the edge. The largest 
differences at the hot end are 14.8 K and 33.3 K for the low and high mode, respectively 
(heater at 23 V in both cases). The differences of the center and edge temperatures in the 
ambient side are approximately 6K for both cases. This spread is represented by 
horizontal error bars shown in Figs. 4.5 and 4.6. The other possible source of the 
discrepancies in resonance frequency and 1/Q at high AT for the high mode is the entire 
temperature distribution of the prime mover. As shown in Figs. 4.3 and 4.4, the major 
change in temperature profile is concentrated at the region adjacent to the hot heat exchanger 
and stack. It is difficult to accurately characterize the entire temperature distribution of the 
prime mover using only 10 thermocouples. This is an important difference between the 
conventional and annular prime movers. In a conventional prime mover, there is typically a 
hot temperature duct and an ambient temperature duct. In an annular prime mover the two 
temperatures must merge along the duct. 

Figures 4.7(a), (b), and (c) show the mode shapes for the low mode at different 
temperature differences. The stack/heat exchanger assembly is located between 
x = 0.768 m and x = 0.792 m, indicated by the dashed line. In this experiment, the 
driver is located 90° from the center of the stack which corresponds to a position of 
0.183 m. The driving frequencies are the resonance frequencies as computed by the 
MATLAB program and the AT’s are calculated using Eqs. 4.1 and 4.2. 

Figure 4.7(a) shows the mode shape for the low mode for AT = 0 K. In this case, the 
stack/heat exchanger assembly acts simply like a constriction and there is a velocity 
antinode near the center of the stack/heat exchanger assembly. The agreement between the 


measured and calculated mode shape is good. 
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Figures 4.7(b) and (c) show the mode shapes for the low mode at temperature 
differences of 80 K and 200 K, respectively. The discrepancy between the measured and 
computed mode shape increases as the temperature difference increases. The discrepancy 
appears as a shift in mode shape. However, in both cases, the measured and predicted 
mode shapes show the same general features. In particular, the standing wave ratio 
decreases with increasing AT. Even with this obvious diminishing in the-standing wave 
ratio, 1/Q only changes by 15% from AT = 0 K to AT = 220 K. 

Figure 4.8(a) show the mode shape for the high mode with AT = 0 K. The agreement 
between the measured and calculated mode shape is excellent. As one might expect, there 
is a pressure antinode near the center of the stack/heat exchanger assembly. Figures 4.8(b) 
and (c) show the mode shapes for high mode at temperature difference of 76 K and 192 K, 
respectively. The agreement between measured and calculated mode shapes are better than 
those of the low mode. In contrast to the low mode shapes, the high mode shapes do not 
change appreciably as AT increases. In particular, the standing wave ratio remains high. 
This difference between the two modes may be caused by the discontinuity in the acoustic 
impedance near the end of the stack. For the low mode, there is a velocity antinode near 
the center of the stack which results in a low acoustic impedance. However, the acoustic 
impedance of the stack is much larger for the high mode. As a result, as AT increases, the 
impedance discontinuity should have more effect on the low mode than the high mode. 

The most important result of this section is that the annular prime mover will not reach 
onset under easily attainable conditions The reason is that neither mode shape supports 
thermoacoustic growth. In other words, the low mode has a velocity antinode centered on 
the stack/heat exchanger assembly whereas the high mode has a pressure antinode centered 


on the stack/heat exchanger assembly. It was also found that the standing wave ratio of the 
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low mode diminishes as AT increases. However, this alteration in mode shape does not 


correlate with increased attenuation in the prime mover. 


78 


Temperature of the prime mover stack for low mode 


300 


T2, center of hot end of stack 

T3, center of ambient end of stack 

T4, edge of ambient end of stack 
250 T5, edge of hot end of stack 


200 


Temperature (C) 
on 
© 


100 


50 





0 5 10 15 20 25 
Heater voltage (V) 


Figure 4.1 Temperature of the prime mover stack vs. the heater voltage (low 
mode). 


a2) 


Temperature of the prime mover stack for high mode 
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Figure 4.2 Temperature of the prime mover stack vs. the heater voltage (high 
mode). 
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Temperature distribution of the low mode of the annular prime mover 
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Figure 4.3 The measured temperature distribution of the prime mover (low 
mode, heater at 23 V). The stack is located to the nght of the dotted line. 
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Temperature distribution of the high mode the annular prime mover 
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Figure 4.4 The measured temperature distribution of the prime mover (high 
mode, heater at 23 V). The stack is located to the nght of the dotted line. 
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Figure 4.5 Comparison of calculated and measured resonance frequencies vs. AT 


for the annular prime mover. 
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1/Q vs DeltaT for the annular prime mover 


0.02 
0.015 
0.01% 
0.005 Oo! eee aoa ee Mi Cie 


—— Calculated High mode 
* Measured High Mode 





0 50 100 150 200 
Delta T (Kelvin) 


Figure 4.6 Comparison of calculated and measured 1/Q vs. AT for the annular 


prime mover. 
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Low mode with dT = 0 Kelvin ; freq = 423.1+5.7351 
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Figure 4.7(a) Mode shape for the low mode of the annular prime mover when the 


driver is located 90° from the stack and AT = 0 K. 
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Low mode with dT = 79.63 Kelvin ; freq = 427.445.7221 
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Figure 4.7(b) Mode shape for the low mode of the annular prime mover when 


the driver is located 90° from the stack and AT = 80 K. 
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Low mode with dT = 199.6 Kelvin ; freq = 438.3+5.442) 
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Figure 4.7(c) Mode shape for the low mode of the annular prime mover when the 


driver is located 90° from the stack and AT = 200 K. 
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High mode with dT = 0 Kelvin ; freq = 438.3+2.099i 


= 
op) 


© 
i 


Normalized pressure 
© 
12) 


© 
a) 


a 
nN 


—— Matlab program 
© Measured data 
0.1 — — stack region 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
Position, m 


Figure 4.8(a) Mode shape for the high mode of the annular prime mover when 


the driver is located 180° from the stack and AT = O K. 
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High mode with dT = 75.73 Kelvin ; freq = 442.74+2.245i 
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Figure 4.8(b) Mode shape for the high mode of the annular prime mover when 
the driver is located 180° from the stack and AT = 76 K. 
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Figure 4.8(c) Mode shape for the high mode of the annular prime mover when 


the driver is located 180° from the stack and AT = 192 K. 
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B. CONSTRICTED ANNULAR PRIME MOVER 


In this section, results for the annular prime mover with a constriction located 90° 
from the center of the stack/heat exchanger assembly are presented. The constriction is 45° 
long and three different porosities 0.7, 0.3, and 0.1 are used. First, the comparisons of the 
measured and calculated resonance frequencies of the constricted’ prime mover are 
presented. Next, the comparisons of the measured and calculated Q of the constricted 
prime mover are discussed. The end corrections are included in the calculations using the 
two methods discussed in Chapter II. This section concludes with the comparisons of the 
measured and computed mode shapes for the constricted prime mover. 

As mentioned earlier, the measured resonance frequencies and Q’s are determined 
from pole-zero analysis of the frequency response. In the measurement, the driver is 
located 45° from the center of the stack assembly and the constriction. The theoretical 
resonance frequencies and Q’s were calculated for four different cases, represented by A, 
B, C, and D. All the results are listed in Tables 4.1 and 4.2. In these tables, Column A 
represents the calculated results without any end corrections. Columns B and C represent 
the results with end corrections for the constriction using the conformal transformation 
method and higher order mode method, respectively. It is worth pointing out that there are 
no appreciable differences between of the results for 1/Q calculated from the conformal 
mapping transformation and the higher order mode method when the end corrections in the 
constriction are included. Column D gives the results for the case of applying end 
corrections to both the constriction and the stack using the higher order mode method. The 
end corrections in the prime mover stack are approximated by applying the end corrections 
to one slit and summing the impedance of the slits in parallel. The temperature differences 


were determined by Eqs. (4.1) and (4.2). Table 4.1 lists the calculated and measured 


| 


resonance frequencies for the low and high modes for the various constrictions. Table 4.2 
lists the calculated and measured Q for the low and high modes for the various 
constrictions. In these tables, Rs represents the area ratio of the constriction, i.e., the 
porosity. 

It is seen that the inclusion of end corrections has significant effects on the predictions 
for the low mode. For example, in the case of Rs =0.1, the averaged differences 
between the calculated and measured resonance frequencies are 1.94 % (no end 
corrections), 0.67% (end corrections using conformal transformation method), and 0.39% 
(end corrections using higher order mode method). On the contrary, including end 
corrections for the constriction has little effect on the predictions for the high mode. This 
point is evident from Figs. 4.9 and 4.10, where it is seen that inclusion of end corrections 
in the constriction is important for the low mode only. However, inclusions of the end 
corrections in both the constriction and stack has more effect on the high mode than on the 
low mode. 

The point to be taken from these figures is that the program does a good job of 
predicting the resonance frequencies, even in the absence of including end corrections. We 
have confidence in the end corrections applied to the constriction. This confidence is the 
result of preliminary work by Choe (Choe, 1997). Application of end corrections to the 
stack is only intended to be approximate, to indicate that they should have an effect. The 
contribution of this work is to point out the necessity of including end corrections. A 
detailed analysis of stack end corrections is an area for future work. 

Figure 4.11 shows the comparisons of the calculated and measured 1/Q of the low 
mode for the unconstricted and constricted prime mover with various porosities. In this set 
of experiments, when the constriction porosity is 0.1, the prime mover reaches onset at a 


temperature difference of 260 K. The horizontal error bar in Fig. 4.11 indicates the 
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uncertainty in determining the onset temperature from the measure temperature distribution 
as discussed previously. In all cases, the calculated results under estimate 1/Q. Therefore, 
the calculated results show a significantly lower onset value of AT in the case of Rs = 0.1. 
The reason for this discrepancy is not fully understood. Nonetheless, both the measured 
and calculated 1/Q indicate the same trends. First, the relative values for 1/Q are consistent. 
Also, both the predicted and measured results for the constricted prime mover intersect 
within a narrow range of AT. The predicted values intersect at AT = 100K and the 
measured data intersect at AT = 175 K. The 1/Q curve for Rs =0.1 begins with the 
highest value of the three constricted cases and ends with the lowest value after the 
intersection at AT =100 K. In contrast, 1/Q for Rs = 0.7 begins with the lowest value 
and ends with the largest value at higher AT. 1/Q for Rs =0.3 remains in between 1/Q 
for the other two porosities. 

The conclusion to be drawn from Fig. 4.11 is that it is predicted that a constricted 
annular prime mover will reach onset, and it does. Although, the agreement between the 
measured and calculated results are only fair, there is very good qualitative agreement. 
This indicates that although there are still some details to be investigated, the general 
behavior is understood. 

Figure 4.12 shows the comparisons between the calculated and measured 1/Q for the 
high mode for the unconstricted and constricted prime mover. It is seen that the 
temperature difference has less effect on 1/Q for the high mode than for the low mode. 
Over the temperature span from AT =0 K to AT = 200 K, the change in 1/Q is less than 
30 % in all cases. The same comments about the general qualitative agreement apply to the 
high mode results as they do to the low mode results. 

The question to be addressed now is why does a constricted annular prime mover 


reach onset. The answer lies in the effects that a constriction has on mode shapes. 
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Figures 4.13(a) and (b) show the comparison of the measured and computed mode 
shapes of the low mode for Rs = 0.1. In this case the constriction is much longer than the 
stack assembly (Lyon crriction!“stack= 3-29) and the porosity of the constriction is much smaller 
than that of the stack assembly. As a result, the mode shape is dominated by the 
constriction. There is a pressure node near the center of the constriction and a pressure 
antinode between the constriction and the stack assembly nearest (but not at) the hot end of 
the stack. End corrections for the constriction are included using the higher order mode 
method. Figure 4.13(a) is for AT=OK _ and Fig. 4.13(b) is for AT =288 K. The 
frequencies are the calculated values. It is observed that there is a good agreement between 
the measured and predicted mode shape at AT =0 K. This mode shape is favorable to 
thermoacoustic growth. The agreement in mode shape worsens at AT =288 K. It should 
be noticed that the imaginary part of the calculated resonance frequency is negative, which 
suggests the onset has been reached. 

Figures 4.14(a) and (b) show the comparison of the measured and computed mode 
shapes of the high mode for Rs =0.1. The agreement between the measured and 
calculated mode shape is good. As one would expect, there is a pressure antinode near the 
center of the constriction. 

One important common feature of Figs. 4.13 and 4.14 is that the mode shapes are 
essentially determined by the constriction and are independent of the temperature 
difference. This is a key to making the annular prime mover reach onset. 

The main conclusion from this section is that to build a functional annular prime 
mover, the relative locations of the constriction and the stack, and their relative porosities 
are the crucial factors. This is equivalent to incorporating some type of dominating 


boundary conditions into the prime mover to alter the mode shapes. 
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Table 4.1 Comparisons of the calculated and measured resonance frequencies for the 
constricted annular prime mover. Case A: without end corrections; Cases B, C, D: 
End corrections with different methods. 
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231 S165 34 34.1 22.8 31.2 
Table 4.2 Comparisons of the calculated and measured Q for the constricted 
annular prime mover. Case A: without end corrections; Cases B, C, D: End 


corrections with different methods. 
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Figure 4.9 Comparisons of the measured and calculated resonance frequencies of 
the low mode for the constricted prime mover with a porosity of 0.1. Method A 


is the conformal mapping transformation and Method B is the higher order mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.1 
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Figure 4.10 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with a porosity of 0.1. Method 


A is the conformal mapping transformation and Method B is the higher order 
mode. 
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Figure 4.11 Comparisons of the measured and predicted 1/Q of the low mode for 


the unconstricted and constricted prime mover. Symbol o represents the 
measured 1/Q of unconstricted prime mover. Symbols *, +, and x represent the 


measured 1/Q for the constricted prime mover with porosities of 0.1, 0.3, and 
0.7, respectively. The lines represent the predicted values. No end corrections 


are included in the predicted values for the unconstricted prime mover. 
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1/Q for the high mode of the annular prime mover 
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Figure 4.12 Comparisons of the measured and predicted 1/Q of the high mode 


for the unconstricted and constricted prime mover. Symbol o represents the 
measured 1/Q of unconstricted prime mover. Symbols *, +, and xX represent the 


measured 1/Q for the constricted prime mover with porosities of 0.1, 0.3, and 
0.7, respectively. The lines represent the predicted values. No end corrections 


are included in the predicted values for the unconstricted prime mover. 
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Figure 4.13(a) Mode shape of the low mode of the constricted prime mover 

(Rs = 0.1) when the driver is located 45° from the stack and AT=0K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Low mode with end corrections ; dT = 228.2 Kelvin ; freq = 307.4-1.165i ; Rs = 0.1 
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Figure 4.13(b) Mode shape of the low mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 228 K. 


The calculated results are based on the higher order mode method. The frequency 
is the calculated value. 
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Figure 4.14(a) Mode shape of the high mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT=0OK. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Figure 4.14(b) Mode shape of the high mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 228 K. 


The calculated results are based on the higher order mode method. The frequency 


is the calculated value. 
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C. TWO STACK ANNULAR PRIME MOVER 


The knowledge gained from the last two sections leads to the design of a symmetric 
two stack annular prime mover in which the constriction in the last section is replaced by a 
stack assembly. This proposed design of the annular prime mover, as shown in Fig. 4.15, 
has two identical stack/heat exchanger assemblies. The two stack assemblies are placed 
such that the hot heat exchangers face each other. 

The MATLAB program was modified to predict the performance of this two stack 
annular prime mover, by adding one more guess of the time-averaged energy flux H 2 
along the second stack to match the temperature of the hot heat exchanger. Two different 
cases are studied here. In case | the duct between the two heat exchangers is held at 
ambient temperature T, = 293 K while in case 2 it is held at the desired temperature of the 
hot heat exchanger 7,. The reason for studying case | is that the current apparatus can be 
easily modified to conduct the experiment. However, because of the temperature 
discontinuity near the junction of the duct and hot heat exchanger, it is difficult to model 
this case theoretically and computationally. It is known that the temperature discontinuity 
in the duct will also create a small discontinuity in the work flow (Rott, 1969). On the 
contrary, although the theoretical analysis of case 2 is simpler than case 1, investigating it 
requires a new apparatus. 

Figure 4.16 shows the calculated resonance frequencies of the two-stack prime mover 
for the two cases. In case 1, it is seen that the resonance frequencies for the high mode 
increase approximately linearly with temperature difference. However, AT has little effect 
on the resonance frequencies of the low mode in case 1. The resonance frequencies for 
both the low and high mode in case 2 increase approximately linearly with temperature 


difference. Because the duct between the two hot heat exchangers is held at 7,, the 
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increase in resonance frequency is bigger than that in case 1. It should be pointed out that 
in case | at very low temperature difference (AT < 50 K), the resonance frequencies of 
the low mode and high mode are very nearly equal. The MATLAB program always 
converges to only one mode. Forexample, at AT = 0 K, the program always converges 
to the high mode, while at AT = 25 K_ it always converges to the low mode. As AT 
increases, the modes are well separated and no such condition was observed. This 
condition is not present in case 2. It was always possible to distinguish between the low 
mode and high mode. It should also be pointed out that for case 1 at AT=0OK, the 
resonance frequency for the high mode is 424.6 Hz and is 424.0 Hz for the low mode. 
The results indicate the possibility of mode level repulsion (Pippard, 1985 and 1989; 
Laraza and Denardo 1997). 

Figure 4.17 shows the calculated 1/Q of the low mode and high mode for the two- 
stack prime mover for the two cases. It is evident that the high mode reaches onset at 
AT = 190 K for case 1. The onset temperature difference for the high mode in case 2 is 
approximately 350 K. An interesting feature different from the 1/Q for the constricted 
prime mover is the linearly increase and decrease in 1/Q for the low mode and high mode, 
respectively. Similar dependence has been observed for a rigid-rigid prime mover (Lin, 
1989; Atchley, et al., 1992). 

Frequency splitting in the two-stack annular prime mover offers the possibility that the 
harmonics of the fundamental frequency and the overtones of the resonator will not 
coincide, thus detuning the system. This means that if high amplitudes are achieved above 
onset, waveform distortion and shocking should not be a problem (Coppens, 1971; 
Atchley and Hofler, 1990; Atchley et. al. 1993). 

Figures 4.18 to 4.21 show the calculated mode shapes for the low mode and high 


mode for the two cases. As one would expect from symmetry, there is a pressure antinode 
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between the two stacks for the high mode and a pressure node between the two stack 
assembly for the low mode. It is evident that only the high modes in both cases support 


thermoacoustic growth, as is confirmed by Fig. 4.17. 
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Figure 4.15 The two-stack annular thermoacoustic prime mover. 
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Figure 4.16 The calculated resonance frequencies for the two stack annular prime 
mover. Case 1: The duct between the two hot heat exchangers is held at T,,. 
Case 2: The duct between the two hot heat exchangers is held at 7,. 
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Calculated 1/Q vs. dT for the two stack annular prime mover 
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Figure 4.17 The calculated 1/Q for the low mode and high mode for the two 
stack annular prime mover. Case 1: The duct between the two hot heat 


exchangers is held at 7,. Case 2: The duct between the two hot heat exchangers 
is held at T,,. 
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Mode shape for of the two stack prime mover, dT = 200K; freq = 426.1+14.74i 
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Figure 4.18 The calculated mode shape of the low mode of the two stack annular 


prime mover at AT = 200 K. Case 1: The duct between the two hot heat 
exchangers is held at 7, = 293 K. 
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Mode shape for of the two stack prime mover, dT = 199.2K; freq = 433.5-0.614i 
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Figure 4.19 The calculated mode shape of the high mode of the two stack annular 
prime mover at AT = 200 K. Case 1: The duct between the two hot heat 
exchangers is held at 7, = 293 K. 
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Figure 4.20 The calculated mode shape of the low mode of the two stack annular 


prime mover at AT = 200 K. Case 2: The duct between the two hot heat 
exchangers is held at T, = 493 K. 
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Figure 4.21 The calculated mode shape of the high mode of the two stack annular 
prime mover at AT = 200 K. Case 2: The duct between the two hot heat 


exchangers is held at T, = 493 K. 
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V. SUMMARY AND CONCLUSION 


The objective of this dissertation is to investigate a thermoacoustic prime mover having 
periodic boundary conditions both experimentally and theoretically. The experimental 
approach has been to design, build, and test an annular thermoacoustic prime mover. The 
resonance frequencies, quality factors, and the mode shapes of the prime mover are 
measured as functions of the temperature difference across the prime mover stack. The 
experimental resonance frequencies and Q’s are determined from pole-zero analysis of the 


measured frequency response. 


An important aspect of this dissertation is the computational analysis of a 
thermoacoustic device with periodic boundary conditions. Two approaches were taken. 
First an existing program, DeltaE, was applied to this novel geometry. It was found that 
DeltaE can be applied to the problem, but has some disadvantages. It is difficult to extract 
O’s below onset and awkward to extract mode shapes. It is also difficult to include 


temperature gradients in ducts. 


Because of these problems, a new program was developed. It uses the same basic 
approach as DeltaE, but overcomes the disadvantages. The complex eigenfrequency is 


used to extract resonance frequencies and Q’s. 


The results of this dissertation indicate that under all but perhaps extreme conditions an 
annular prime mover will not reach onset. The reason is well understood. The eigenmodes 
of a single-stack annular prime mover do not support thermoacoustic growth. This finding 
lead to the investigation of a constricted annular prime mover. The idea is that the 


constriction could impose dominating boundary conditions on the prime mover thus 
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altering the eigenmodes. It was shown that decreasing the porosity of the constriction 
forced the prime mover into onset. The porosity and location of the constriction is used to 


adjust the eigenmodes of the system so that thermoacoustic growth is supported. 


Another way to produce eigenmodes that support thermoacoustic growth is to 
introduce a symmetry into the system which shifts the nodes and antinodes of the modes 
off the stacks, such that for one of the modes a pressure antinode is _ the hot end of the 
stacks. It is predicted that a two-stack annular prime mover will reach onset. This design 


is an area for future research. 


Another area for future study is end correction for stacks. It was shown in this 
dissertation that inclusion of end corrections does affect the predicted resonance frequencies 
and Q’s. Determining the proper end corrections for realistic stack geometries is likely to 


be a nontrivial project. 
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APPENDIX A.1 THE DELTAE INPUT FILE FOR THE LOW FREQUENCY 


MODE OF THE ANNULAR THERMOACOUSTIC PRIME MOVER 


TITLE Annular resonator (file : to2dTflowC.in) 
!Created @23:18:35 02-MAY-97 with DeltaE Vers. 3.0b1 for the Power Macintosh 


|_---..--------------------------- Q ----------------------- 

BEGIN Initial 

1.0130E+05 aMeanP Pa 16.63 

427.0 b Freq. Hz -43.90 

293.0 cT-beg K 6.7062E-04 

16.63 dlpl@0 = Pa G 38.50 

-43.90 e Ph(p)O deg G_ -1.8309E-03 

6.7062E-04 fIUI@O0 m3/s G 

38.50 g Ph(U)O deg G 

air Gas type 

ideal Solid type 

|--- ------------------------------ | ----------------------- 

SOFTEND 1 

0.0000 a Re(Z) (t) 16.63 

0.0000 b Im(Z) (t)  -43.90 
6.7062E-04 
38.50 
7.3673E-04 
7.3673E-04 
2.0847E-02 

sameas 0 Gas type -0.1564 

ideal Solid type Zou 

|--------------------------------- 2 -n---2---2-------------- 

ISODUCT Duct 

2.6300E-03 aArea m2 42.05 

0.2050 bPeim m -48.73 

3.1788E-02 cLength m 6.2404E-04 
38.20 

sameas 0 Gas type 7.0088E-04 
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A |lpl@0 G( 0d) P 
B Ph(p)0 G( Oe) P 
CIUI@O G(Of) P 
D Ph(U)0 G( 0g) P 
E HeatIn G(27e) P 


A Ipl Pa 
BPh(p) deg 
ClUl = m43/s 
DPh(U) deg 
E Hdot W 
FWork W 
G Re(Z) 

H Im(Z) 

I T K 

A Ipl Pa 
BPh(p) deg 
CIUl = m43/s 
D Ph(U) deg 
EHdot W 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas QO 


ideal 


{SODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


Duct 
a Area 
b Perim 


c Length 
Gas type 


Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


c Length 


Gas type 
Solid type 


c Length 


sameas 0 Gas type 


ideal 


[ISODUCT 
2.6300E-03 


Solid type 


Duct 


a Area 


m 


mm 


m2 


7.0088E-04 


64.95 
-49.96 

5.3889E-04 
37.86 
6.6858E-04 
6.6858E-04 


83.85 
-50.55 

4.2045E-04 
3736 
6.4101E-04 
6.4101E-04 


97.57 
-50.93 

2.7606E-04 
36.44 
6.1820E-04 
6.1820E-04 


105.3 
-51.21 
1.1484E-04 
33.10 
5.9903E-04 
5.9903E-04 
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F Work 


A Ipl 

B Ph(p) 
C1IUI 

D Ph(U) 
E Hdot 
F Work 


A Ipl 

B Ph(p) 
C {Ul 

D Ph(U) 
E Hdot 
F Work 


A Ipl 

B Ph(p) 
C IU 

D Ph(U) 
E Hdot 
F Work 


A Ipl 

B Ph(p) 
C IU 

D Ph(U) 
E Hdot 
F Work 


W 


Pa 
deg 
m‘3/s 
deg 
W 
W 


Pa 
deg 
m‘3/s 
deg 
Ww 
Ww 


Pa 
deg 
m‘3/s 
deg 
W 
W 


Pa 
deg 
m‘3/s 
deg 
WwW 
WwW 


0.2050 
3.1788E-02 


sameas Q 


ideal 


b Perim m 


c Length m 


Gas type 
Solid type 


ISODUCT Duct 


2.6300E-03 aArea m2 
0.2050 bPerim m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
eee - 2 __-.___-_ 
ISODUCT Duct v 
2.6300E-03 aArea m2 
0.2050 bPerim =m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
i 
ISODUCT Duct 8 
2.6300E-03 aArea m2 
0.2050 bPerim =m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 


ISODUCT Duct 


2.6300E-03 
0.2050 
3.1788E-02 


sameas QO 


a Area m2 
b Perim m 


c Length m 


Gas type 


-51.45 
5.5371E-05 
-130.1 

5.8148E-04 

5.8148E-04 


101.0 
-51.67 

2.1977E-04 
-138.8 
5.6316E-04 
5.6316E-04 


89.39 
-51.92 

3.7131E-04 
-140.0 
5.4183E-04 
5.4183E-04 


1222 
-52.23 

4.9997E-04 
-140.6 
5.1603E-04 
5.1603E-04 


309 
-52.75 

5.9773E-04 
-140.9 
4.8535E-04 
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B Ph(p) deg 
CIU! m‘3/s 
DPh(U) deg 
EHdot W 

F Work W 


A Ipl Pa 

B Ph(p) deg 
CIlUl = m/43/s 

DPh(U) deg 
EHdot W 

F Work W 


A Ipl Pa 
BPh(p) deg 
CIUlL = m43/s 
DPh(U) deg 
EHdot W 
FWork W 


A Ip! Pa 
BPh(p) deg 
ClU! = m43/s 
D Ph(U) deg 
EHdot W 
FWork W 


A Ipl Pa 
BPh(p) deg 
CIUl = m43/s 
DPh(U) deg 
EHdot W 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 


Duct 
a Area 
b Perim 


c Length 


Gas type 
Solid type 


b Perim 


c Length 


Gas type 
Solid type 


c Length 


Gas type 
Solid type 


c Length 


Gas type 
Solid type 


m 


m 


4.8535E-04 


25.84 
-54.17 

6.5855E-04 
-141.1 
4.5059E-04 
4.5059E-04 


1.340 
-165.9 

6.7865E-04 
-141.3 
4.1353E-04 
4.1353E-04 


26.91 

131.0 

6.5679E-04 
-141.5 
3.7652E-04 
3.7652E-04 


51256 

129.7 

5.9432E-04 
-141.6 
3.4190E-04 
3.4190E-04 
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eg e @ 2 e@ 2 2 SO eS Oe OSB oO ee BS eS eS SC Ses eee TSF TE 


FWork W 
A Ipl Pa 
BPh(p) deg 
ClUl = m/43/s 
DPh(U) deg 
EHdot W 
FWork W 

A Ip! Pa 
BPh(p) deg 
CU! = m3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUl = m/43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
Crm me3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 


0.2050 


b Perim m 


3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
|--------------------------------- 17 
ISODUCT Duct 

2.6300E-03 aArea m2 
0.2050 bPerim m 
3.1788E-02 cLength m 
Sameas 0 Gas type 

ideal Solid type 
|--------------------------------- 18 
ISODUCT Duct 

2.6300E-03 aArea m2 
0.2050 b Perim m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 

| --------------------------------- 19 
ISODUCT Duct 

2.6300E-03 aArea m2 
0.2050 b Perim = m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
|--------------------------------- 20 


VDUCER Driver 


1 .QOOOE-09 
0.0000 

1 .QOOOE+04 
0.0000 
-1.0000E+04 


a Re(Ze) ohms 

b Im(Ze) 

c Re(T1) V-s/m‘3 
d Im(T1) V-s/m43 
e Re(T2) Pa/A 


12932 

4.9511E-04 
-141.8 
3.1142E-04 
3.1142E-04 


S899 

128.9 

3.6528E-04 
-142.1 
2.8383E-04 
2.8583E-04 


101.4 
128.8 
2.1289E-04 

-142.6 
2.6468E-04 
2.6468E-04 


106.5 
128.7 
4.7479E-05 

-146.9 
2.4644E-04 
2.4644E-04 


106.5 

128.7 
1.4217E-04 

-169.5 
3.5727E-03 
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BPh(p) deg 
ClU! = =m43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 

B Ph(p) deg 
ClUl = m43/s 
DPh(U) deg 
EHdot W 
FWork W 

A Ipl Pa 
BPh(p) deg 
ClU! = m43/s 
DPh(U) deg 
EHdot W 

F Work W 

A Ipl Pa 
BPh(p) deg 
CIV] = m/43/s 
DPh(U) deg 
EHdot W 
FWork W 

A Ip! Page 
BPh(p) deg 
CjUl ame3/s 
DPh(U) deg 
EHdot W 


0.0000 
1.QOO0OE-09 
1.QOOOE-09 
1.000 


sameas 0 


ideal 


fIm(T2) Pa/A 


g Re(Zm) Pa-s/m‘3 
h Im(Zm) Pa-s/m‘3 


i AplVol  V 


Gas type 
Solid type 


ISODUCT Duct 


2.6300E-03 
0.2050 
3.1788E-02 


sameas O 


ideal 


a Area m2 


b Perim m 


cLength m 


Gas type 
Solid type 


SODUCT Duct 


2.6300E-03 
0.2050 
3.1788E-02 


sameas QO 


ideal 


2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


2.6300E-03 
0.2050 
3.1788E-02 


ISODUCT Duct 


ISODUCT Duct 


a Area m2 


b Perim m 


c Length m 


Gas type 
Solid type 


a Area m2 


b Perim m 


c Length m 


Gas type 
Solid type 


a Area a4 
b Perim m 


c Length m 


3.5727E-03 
3.3263E-03 
1.000 
1.065 1E-02 
-51.35 
1.QOO0E-04 
-180.0 


108.1 
| 
7.9112E-05 
93.49 
3.5543E-03 
3.5543E-03 


103.1 
125.9 
2.2143E-04 
53.91 
3.5352E-03 
3.5352E-03 


91.74 

124.2 

3.7161E-04 
46.13 
3.5132E-03 
3.5132E-03 


74.83 
122.0 
5.0174E-04 
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FWork W 
G WorkIn W 
H Volts V 

I Amps V 

J Ph(Ze) deg 
K lUxl = m‘3/s 
L Ph(-Ux deg 
A Ipl Pa 
BPh(p) deg 
Gil 8 me3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUl = m‘3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUl = m/43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
CIUIl = m3/s 


sameas 0 


ideal 


ISODUCT Duct 


2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT Duct 


2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


HXFRST Ambient HX 


sameas 3a 
0.6290 
5.0800E-03 
1.7018E-03 


-1.8309E-03 


293.0 
sameas QO 


copper 


STKSLAB Prime Mover Stack 


2.6300E-03 
0.6220 

2.4000E-02 
3.0480E-04 
7.5000E-05 


42.72 
3.4867E-03 
3.4867E-03 


53.48 

118.2 

6.0191E-04 
40.63 
3.4553E-03 
3.4553E-03 


29.40 

108.6 

6.6543E-04 
39.08 
3.4195E-03 
3.4195E-03 


23:02 
103.2 
6.6930E-04 
38.94 
1.5886E-03 

3.3455E-03 
-1.8309E-03 
293.0 


16.26 
-43.67 
6.7080E-04 
38.51 
1.5886E-03 
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DPh(U) deg 
EHdot W 
FWork W 

A Ipl Pa 
BPh(p) deg 
CIUl = m%3/s 
DPh(U) deg 
EHdot W 
FWork W 

A Ipl Pa 
BPh(p) deg 
ClIUl = m%3/s 
D Ph(U) deg 
EHdot W 
FWork W 
A Ipl Pa 

B Ph(p) deg 
CIUl = m/43/s 
DPh(U) deg 
EHdot W 

F Work W 
G Heat W 
H MetalT K 
A Ipl ma 

B Ph(p) deg 
CIUl m%3/s 
D Ph(U) = deg 
EHdot W 


sameas OQ 


kapton 


HXLAST Cold HX 


sameas 3a 
0.7050 
3.0480E-04 
1.2540E-03 
0.0000 
293.0 
sameas 0 


copper 


sameas 0 


ideal 


0.0000 
1A b DlAddr 
30A c D2Addr 


DIFFTAR 
0.0000 

1B b DlAddr 
30B c D2Addr 


Ph(p) mismatch 


7.4220E-04 

293.0 

295.0 
-2.6033E-03 


16.63 
-43.90 

6.7062E-04 
335 1 
7.3668E-04 
7.3668E-04 


-8.5195E-04 


2980 


16.63 


-43.90 


6.7062E-04 
38.51 

7.3668E-04 
7.3668E-04 
2.0846E-02 


-0.1564 


293.0 


8.0065E-05 


=32A? 3.3780E-04 
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ef 2 28 eS SO SC SO SB SSE SB OS OS SO See SE eS Se SS Se ee oS eS ee 


F Work W 
GT-beg K 
H T-end K 
TStkWrk W 
A Ipl Pa 
BPh(p) deg 
Clu! m3/s 
DPh(U) deg 
EHdot W 
F Work W 
GHeat W 
H MetalT K 
A Ipl Pa 
BPh(p) deg 
ClUlL m43/s 
DPh(U) deg 
EHdot W 
F Work W 
G Re(Z) 

H Im(Z) 

ite |; K 

A D1-D2 

A D1-D2 


DIFFTAR [Ul mismatch 


0.0000 a TargDi =33A? 1.3884E-09 A D1-D2 
1C b DlAddr 

30C c D2Addr 

|--------------------------------- 34 --------------------------------- 
DIFFTAR Ph(U) mismatch 

0.0000 a TargD1 =34A? -1.0691E-04 A D1-D2 


1D b D1Addr 
30D c D2Addr 


! The restart information below was generated by a previous run 
! You may wish to delete this information before starting a run 
! where you will (interactively) specify a different iteration 
! mode. Edit this table only if you really know your model! 
INVARS 50405 0 60 727 5 
TARGS 529 631 132 133 134 1 
SPECIALS 0 
PLTVAR 702-110102030405 7120 1 
! Plot start, end, and step values. May be edited if you wish. 
! Outer Loop: | Inner Loop . 
410.0 427.0 0.2000 1.000 1.000 1.000 
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APPENDIX A.2 THE DELTAE INPUT FILE FOR THE HIGH FREQUENCY MODE 
OF THE ANNULAR THERMOACOUSTIC PRIME MOVER 


TITLE Annular resonator (file :: to2dTfhighC.in) 
!Created @ 12:00:08 02-MAY-97 with DeltaE Vers. 3.0b1 for the Power Macintosh 
|--------------------------------- Q) --------------------------------- 
BEGIN Initial 
1.0130E+05 aMeanP Pa 0.0000 A lpi@O G( 0d) P 
430.5 b Freq. Hz 0.0000  B Ph(p)0 G( 0e) P 
293.0 cT-beg K 0.0000 CIUI@0O Gi Of) P 
500.0 dipi@O Pa G 0.0000 DPh(U)0 G( Og) P 
0.0000 ePh(p)0 deg G 0.0000 EHeatIn G(27e) P 
2.9000E-04 fIVI@O m3/s G 
0.0000 g Ph(U)O deg G 
air Gas type 
ideal Solid type 
i ae ee a es 
SOFTEND 1 
0.0000 a Re(Z) (t) 0.0000 =A Ipi Pa 
0.0000 b Im(Z) (t) 0.0000 BPh(p)_ deg 
0.0000 CIiUl m3/s 
0.0000 DPh(U)~ deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 G Re(Z) 
sameas Q Gas type 0.0000 4H Im(Z) 
ideal Solid type 0.0000 I T K 
ee ere ae ee en eee hemes ee nese 
ISODUCT Duct 
2.6300E-03 aArea m2 0.0000 A Ipl Pa 
0.2050 b Perm =m 0.0000 BPh(p)_ deg 
3.1788E-02 cLength m 0.0000 CIUI m3/s 
0.0000 DPh(U)~ deg 
sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
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ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0Q 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 


Duct 
a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 
Gas type 


Solid type 


a Area 


b Perim 


eae eee eee ee eee es = 


0.0000 
0.0000 


A Ipl Pa 
BPh(p) deg 
ClIUl = m43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
ClIUl = =m‘3/s 
DPh(U) deg 
EHdot W 
F Work W 
A Ipl Pa 
BPh(p) deg 
ClUl = m43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
CIUl = m‘%3/s 
D Ph(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
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3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


c Length 


Gas type 
Solid type 


b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


m 


m2 
m 


m 


m42 
m 


m 


0.0000 CIUIL m/‘3/s 
0.0000 DPh(U) deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 Alpi Pa 
0.0000 BPh(p) deg 
0.0000 CIUI m‘3/s 
0.0000 DPh(U)~ deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 A lIpl Pa 
0.0000 BPh(p)_ deg 
0.0000 CIUIL m‘3/s 
0.0000 DPh(U)~ deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 A Ipi Pa 
0.0000 BPh(p)_ deg 
0.0000 CIUIl m43/s 
0.0000 DPh(U) deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 Alpi Pa 
0.0000 BPh(p)_ deg 
0.0000 CIUI m‘3/s 
0.0000 DPh(U)~ deg 
0.0000 EHdot W 
0.0000 FWork W 
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ISODUCT Duct 


2.6300E-03 aArea m/‘2 0.0000 A 'Ipl Pa 
0.2050 bPeim m 0.0000 BPh(p)_ deg 
3.1788E-02 cLength m 0.0000 CIUI m3/s 
0.0000 DPh(U) deg 
sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
|--------------------------------- 13 --------------------------------- 
ISODUCT Duct 
2.6300E-03 aArea m2 0.0000 Alpi Pa 
0.2050 bPerim m 0.0000 BPh(p)_ deg 
3.1788E-02 cLength m 0.0000 CIUI m%3/s 
0.0000 DPh(U) = deg 
sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
|--------------------------------- 14 --------------------------------- 
VDUCER Driver 
1.OO000E-09 aRe(Ze) ohms 0.0000 Alpi Pa 
0.0000 b Im(Ze) 0.0000 BPh(p)_ deg 
1.0000E+04 = c Re(T1) V-s/m‘43 0.0000 CIUI m‘43/s 
0.0000 d Im(T1) V-s/m‘3 0.0000 DPh(U)~ deg 
-1.0000E+04 eRe(T2) Pa/A 0.0000 EHdot W 
0.0000 fIm(T2) Pa/A 0.0000 FWork W 


1.0000E-09 = g Re(Zm) Pa-s/m43 0.0000 GWorkIn W 
1.0000E-09 hIm(Zm) Pa-s/m‘3 0.0000 HVolts V 


1.000 iAplVol V 0.0000 IAmps'  V 
0.0000 JPh(Ze) deg 
sameas 0 Gas type 0.0000 KIUxl m‘43/s 
ideal Solid type 0.0000 LPh(-Ux deg 
|--------------------------------- 15 --------------------------------- 
ISODUCT = Constriction 
2.6300E-03 aArea m2 0.0000 A lIpl Pa 
0.2050 bPeim om 0.0000 BPh(p)_ deg 
3.1788E-02 cLength m 0.0000 CIUI m3/s 


0.0000 DPh(U) deg 
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sameas 0 Gas type 

ideal Solid type 

| --------------------------------- 16 
ISODUCT Duct 7 
2.6300E-03 aArea m2 
0.2050 bPeim om 
3.1788E-02 cLength m 
sameas Q Gas type 

ideal Solid type 

| --------------------------------- 17 
ISODUCT Duct 8 
2.6300E-03 aArea m2 
0.2050 b Perm om 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
|--------------------------------- 18 


ISODUCT Duct 


2.6300E-03 aArea m2 
0.2050 b Perm m 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
|------~-------------------------- 19 


ISODUCT Duct 


2.6300E-03 aArea m2 
0.2050 b Perm mm 
3.1788E-02 cLength m 
sameas 0 Gas type 

ideal Solid type 
|--------------------------------- 20 


ISODUCT  Constriction 


—S—SsQ* 2 Se eee Se SP ee 2 Pee es eS Tee SE Se Se eee ee = 


EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUlL = m/43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
CIUI = m3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUI = m/3/s 
DPh(U) deg 
EHdot W 

F Work W 
A Ipl Pa 
BPh(p) deg 
CIUl = m‘43/s 
DPh(U) deg 
EHdot W 
FWork W 
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2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 
ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas 0 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas OQ 


ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


sameas QO 
ideal 


ISODUCT 
2.6300E-03 
0.2050 
3.1788E-02 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


Gas type 
Solid type 


b Perim 


c Length 


Gas type 
Solid type 


c Length 


Gas type 
Solid type 


a Area 
b Perim 


c Length 


m2 
m 


m 


ee See eee eee Ses oe @ 


0.0000 
0.0000 
0.0000 
0.0000 


A Ipl Pa 
BPh(p) deg 
ClIUl = m/‘3/s 
DPh(U) = deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
ClUl = m‘3/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
ClIUl = m*%3/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
B Ph(p) deg 
ClIUl = m43/s 
DPh(U) deg 
EHdot W 
FWork W 
A Ipl Pa 
BPh(p) deg 
ClIUl = m43/s 
DPh(U) deg 
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sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
|--------------------------------- 25 --------------------------------- 
ISODUCT Duct 
2.6300E-03 aArea m2 0.0000 A ipl Pa 
0.2050 b Perm =m 0.0000 BPh(p) deg 
3.1788E-02 cLength m 0.0000 CIUl m‘3/s 
0.0000 DPh(U)- deg 
sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
|--------------------------------- 26 --------------------------------- 
ISODUCT Duct 
2.6300E-03 aArea m2 0.0000 A 'Ipl Pa 
0.2050 b Perim m 0.0000 BPh(p)_ deg 
3.1788E-02 cLength m 0.0000 CIUl m3/s 
0.0000 DPh(U)~ deg 
sameas 0 Gas type 0.0000 EHdot W 
ideal Solid type 0.0000 FWork W 
|-----..--.----------------------- 2] --------------------------------- 
HXFRST Ambient HX 
sameas 3a a Area m2 23.02 A Ipl Pa 
0.6290 b GasA/A 103.2 BPh(p) deg 
5.0800E-03 cLength m 6.6930E-04 CIUl m‘3/s 
1.7018E-03 dy0 m 38.94 DPh(U) deg 
-1.8309E-03 eHeatIh W G 1.5886E-03 EHdot W 
293.0 fEst-T K_ (t) 3.3455E-03 FWork W 
sameas 0 Gas type -1.8309E-03 GHeat W 
copper Solid type 293:0 H MetalT K 
|--------------------------------- 28 --------------------------------- 
STKSLAB Prime Mover Stack 
2.6300E-03 aArea m2 16.26 A Ipl Pa 
0.6220 b GasA/A -43.67 BPh(p) deg 
2.4000E-02 cLength m 6.7080E-04 CIUl m3/s 
3.0480E-04 dy0 m 38.51 DPh(U) deg 
7.5000E-05 eLplate m 1.5886E-03 EHdot W 
7.4220E-04 FWork W 
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293.0 GT-beg K 


sameas 0 Gas type 29350) H T-end K 
kapton Solid type -2.6033 E-03 1StkWrk W 
|--------------------------------- 29 --------------------------------- 
HXLAST Cold HX 

sameas 3a aArea m/2 16.63 A Ipl Pa 
0.7050 b GasA/A -43.90 BPh(p) deg 
3.0480E-04 cLength m 6.7062E-04 CIUl m3/s 
1.2540E-03 dy0d m 38.51 DPh(U) deg 
0.0000 eHeatIn W_  (t) 7.3668E-04 EHdot W 
293.0 fEst-T K =29H? 7.3668E-04 FWork W 
Ssameas 0 Gas type -8.5195E-04 GHeat W 
copper Solid type 293.0 H MetalT K 
epeemeeeens-.._____---.-----2-=2 24) eee nae 
SOFTEND 

0.0000 a Re(Z) (t) 0.0000 A Ipl Pa 

0.0000 b Im(Z) (t) 0.0000 BPh(p)_ deg 


0.0000 CIUI m3/s 
0.0000 DPh(U)~ deg 
0.0000 EHdot W 
0.0000 FWork W 
0.0000 G Re(Z) 


sameas 0 Gas type 0.0000 H Im(Z) 

ideal Solid type 0.0000 I T K 
|--------------------------------- 3] --------------------------------- 
DIFFTAR Ilp| mismatch 

0.0000 a TargDi =31A? 8.0065E-05 A D1-D2 
LA b DlAddr 

30A c D2Addr 

|--------------------------------- 32 --------------------------------- 
DIFFTAR Ph(p) mismatch 

0.0000 a TargDi =32A? 3.3780E-04 A D1-D2 
1B b Dl Addr 

30B c D2Addr 

| oe BRU daa: Veen eee Saeew reer eee 


DIFFTAR Ul mismatch 
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0.0000 a TargDi =33A? 1.3884E-09 A D1-D2 
IC b D1 Addr 
30C ¢ D2Addr 


DIFFTAR Ph(U) mismatch 

0.0000 a TargDi1 =34A? -1.0691E-04 A D1-D2 
1D b D1Addr 

30D c D2Addr 


! The restart information below was generated by a previous mun 
! You may wish to delete this information before starting a run 

! where you will (interactively) specify a different iteration 

! mode. Edit this table only if you really know your model! 
INVARS 50405 0 6 0 727 5 

TARGS 5 29 631 132 133 134 1 

SPECIALS 0 
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APPENDIX B. THE MATHEMATICA PROGRAM FOR THE LEAST 
SQUARE FIT TO THE FREQUENCY RESPONSE 


(* This mathematica program read the plotting data file generated by *) 

(* DeltaE and perform a least square fit to find the resonance frequency 7) 

ClearAll["Global**"] 

(* Read the data file generated by the DeltaE program. *) 

(* In this example, the data file is flowDT100.de 7) 

datfile=ReadList["flowDT100.de" ,Table[Number, {8}]];(* There are eight columns in *) 
(* the data file *) 

(* Select the data range used in the least square fit, column 1 is the frequencies *) 

(* and column 7 is the pressure amplitude *) 

data=datfile[[Range[2,19],{1,7}]] (* Determine the bandwidth of fitting =) 

(* Initial guesses of A, f,, and Q *) 

soln = {A->150,f,->437,Q->100}; 


Attributes[y]={Listable }; 


y[f_,A_,f,_ .Q_] :=A (f/f,)42 Abs[Cot[Pi (f/f,-I/(2 Q))]/(Pi (f/£,-1/(2 Q)))I; 
Attributes[yy] = { Listable}; 

yy(f_] = y[f,A,f,,Q] /.soln; 

xlist = Transpose[data][f1 ]]; 

ylist = Transpose[data]|[[2]]; 

rmserror[A_,f,_,Q_] := Sqrt[ Apply[Plus,(ylist - y[xlist, A,f,,Q])*2]]; 
rmserror[A,f,,Q]/.soln; 

FindMinimum[rmserror[A,f,,Q],{A, { 100,200} }, { £,, {435,439} },{Q, {90,110} }] 
(* A is an arbitrary constant, f, the resonance frequency, Q the quality factor *) 

(* End of program *) 


S35 
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APPENDIX C. THE MATLAB PROGRAM FOR THE ANNULAR 
THERMOACOUSTIC PRIME MOVER 


Function updateF.m: The function performs the actual iteration to find the solutions for 


the constricted annular prime mover. 


% Begin of function % 
function [y,iterations]=updateF(y); 


% 

% usage [newy,iterations] = updateF(y) 

% This function implements picards method for updating y 
J solving f(y,x)=0 for an annular resonator 


Jo The input y vector is a 8 by 1 vector: [Tm;pre;pim;Ure;Uim;fre;fim;H2dot] 
% only the last 5 elements (i.e. Ure, Uim, fre, fim and H2) are updated 

Jo newy = y - inv(J)*f(y,x); 

% where J is the Jacobian of f(y) 

%o REQUIRES: fdjacF.m, funcvF.m, computeL.m, 

% Note: Before run this program, make sure to adjust the area ratio and select the 
% right temperature (Just active the line for the temperature profile) 

% Declare the global variables 

global T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 Th Tc Tmatrix Rs La Las R Rst r 

r = 2.565e-2; % Equivalent radius of the annular resonator, m 

fs = 0./: % The area ratio for the constriction 


% Compute the equivalent inductance of constriction 


d = 0.0529; % Width of annulus, m 

h = 0.0497; % Height of annulus, m 

dc = d*sqrt(Rs); % Width of the constriction, m 

he = h*sqrt(Rs); % Height of annulus, m 

dst = 4.28e-2; % Width of a slit in the stack m 

hst = 6.12e-4; % Height of a slit in the stack m 

La = computeL(d,h,dc,hc); % Compute the equivalent inductance of constriction 


Las = computeL(d,h,dst,hst); % Compute the equivalent inductance of a slit in the 
J stack 
R = sqrt((dc*hc)./(d.*h)); 
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Rst = sqrt((dst*hst)./(d.*h)); | % Compute the equivalent area ratio for one slit 
% The measured temperature profiles for all the measurements 
% When run this main program, select the profile based on the porosity 
% The low mode and high mode share the same temperature profile 
% These temperatures are in °C 
% Low mode and high mode with constriction Rs=0.1 % Heater voltage 
% Tmatrix=[21.0 21 21 21 21 21 21 21 21 20.6); Coie ater aan) V. 
Vomlbimatiix—(06.9 23.4 22.5 58.8 22.2 21-9 22.2 32.5 25.7 21.4]; % heater = 8V 
% Tmatrix=[100.7 25.8 24.3 86.8 23.8 23.2 23.7 41.2 30.4 22.6]; % heater = 11V 
% Tmatrix=[137.8 28.5 26.2 119.2 25.4 24.4 25.2 51.7 36.2 23.6]; % heater = 14V 
% Tmatrix=[218.6 36.6 32.4 194.5 30.6 28.7 30.2 82.1 51.8 27.3]; % heater = 20V 
% Tmatrix=[278.6 43.9 38 253.1 35.2 32.3 34.5 113.1 66.2 30.2]; % heater = 25V 
% Low Mode and high mode with constriction Rs=0.3 
% Tmatrix=[29.7 29.5 29.4 29.6 29.5 29.4 29.4 29.3 29.5 29.5]; % heater = OV 
% Tmatrix=[74.3 30 29.1 66.8 28.9 28.5 28.8 37.8 32.2 28.4]; % heater = 8V 
% Tmatrix=[106.6 31.1 29.6 93.8 29.1 28.4 29.0 44.5 35.1 28.2]; % heater = 11V, 
% Tmatrix=[143.9 33.1 30.7 126 29.9 28.9 29.7 54.2 39.4 28.4]; % heater = 14V, 
% Tmatrix=[220.7 38.5 34.4 196.2 32.6 30.6 32.1 81.4 50.8 29.5]; % heater = 20V 
% Tmatrix=[278.4 44.9 39.1 251.5 36.3 33.2 35.5 110.3 63.1 31.5]; % heater = 25V 
% Low and high Mode with constriction Rs=0.7 
Comlimathiie-—2 5952 36) 250) 29-5 25-9 25.5 23.5 23.5 23.4); % heater = OV 
% Tmatrix=(74.9 25.3 24.3 65.7 24 23.6 23.9 55.1 23.7 23.4]; % heater = 8V 
% Tmatrix=[108.5 27.2 25.5 94 24.9 24.2 24.8 74.5 30.4 23.9]; % heater = 11V 
% Tmatrix=[147.7 29.8 27.3 128.2 26.3 25.1 26.0 100.1 34.5 24.6]; % heater = 14V 
% Tmatrix=[221.9 35.6 31.2 196 29.4 27 28.8 157.1 43.6 26.2]; % heater = 20V 
Vow Winatiix=(280.7 42.5 36.2 252.9°33.3 29757 32-4e 2 53.5 29:5]; % heater 25 V 
Tmatrix = Tmatrix+273; % Change temperature to Kelvin 
T2=Tmatrix(1); 
T3=Tmatrix(2); 
T4=Tmatrix(3); 
T5=Tmatrix(4); 
T6=Tmatrix(5); 
T7=Tmatrix(6); 
T8=Tmatrix(7); 
T9=Tmatrix(8); 
T10=Tmatrix(9); 
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T11=Tmatrix(10); 
Th = (2*T2+TS)/3; 
Te = (2*T3+T4)/3; 


Nmax = 100; 
y=y(:); 
y(1,1) = Th; 


iterations=0; 


% Weighted hot end targeted temperature, K 

% Weighted cold end temperature, K 

% Maximum number of integration 

% Make sure the input vector y is a column vector 


% Fixed the starting temperature at T,, 


% Main loop for iterations 


while(1) 
iterations = iterations + 1; 
fvec=funcvF(y); % Compute the difference between the targets and % 
% calculated results %o 


checky = norm(fvec./y(1:5)) % Check the norm of the ratio of the difference % 
if checky<1.0E-3 % Set the converging threshold %o 
fprintf(1,‘Done\n’); 


y=y 


% Print out the solution % 


[fvec,yend] = funcvF(y); 


yend = yend' % print out the integration results at the ends % 
Q = yend(1,6)/(2*yend(1,7)) % Find the quality factor % 


elseif iterations == Nmax % Terminate the iteration if iteration = N,,, % 


fprintf(1,'Exceed max # of iterations’); 


return 


else 


i aiajack(y,fvec); % calculate the Jacobian % 


diffy = -J\fvec; 


% Use the LU factorization % 


if checky<1.0e-1 
diffy = diffy/2; % Close to converge, adjust the step size % 


elseif checky<5.0e-2 
diffy = diffy/4; 

elseif checky<1.0e-2 
diffy = diffy/6; 


end 


newy=y 


newy(4:8) = y(4:8)+diffy; % Perform a Newton iteration to update y % 


y=newy; 


% only update Ure, Uim, fre, fim and H2 % 


fprintf(1,'Computing Update, iterations %4.0f \n',iterations) 
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end 
end 
% 
% End of updateF.m % 
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Function funcvF.m: The function integrate the coupled ordinary differential equations of 
the annular prime mover. 


% Begin of function % 

function [fvec,yend] = funcvF(yin) 

Jo 

% Includes the ends effect from high order mode method 

% This function integrates the annular resonator with the stack/heat exchanger assembly 
% at the end of the resonator (including ambient and hot heat exchanger) 
% constriction locates between xcon1 to xcon2 

% Requires: ductfunc2.m, ductfunc21.m, ductcon.m, 

% , stackfun.m, ode45.m, ambhxfunc.m, hothxfunc, zstuff.m 

% Declare the global variables 

global T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 Th Tc Tmatrix Rs La Las R Rst 
yin=yin(:); 

d= 0.0529; % Width of annulus m 

h = 0.0497; % Height of annulus m 


dc = d*sqrt(Rs); % Width of the constriction m 
he = h*sqrt(Rs); % Height of constriction m 
cp=1005; % Isobaric heat capacity per unit mass 


% These are the positions of the thermocouples, m 


xp0= 0; 

xp11=0.00293; % Thermocouple # 9 

xp12=0.084345; % Thermocouple # 10 

xp13=0.39615; % Thermocouple # 11 

xp14=0.707955; % Thermocouple # 7 

xp15=0.7556; % Thermocouple # 8 

mpi =0°76292; % Ambient heat exchanger 

xp2="0:768; % Prime mover stack 

mee = 0.7920; % Hot heat exchanger 

xcon1=0.133863; % 45° long constriction, 90° from the center of the stack assembly 
xcon2= 0.23290; 

xpL=0.7923; % Effective length of the annular resonator, m % 


% These are not the resolution of the integration 


% These are the points in which the results are extracted 
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xn le=320: xn2e=i30: 

Xi = 50. xn4 = 10; 

fvec = zeros(5,1); % Set up the fvec vector 
xllspan = linspace(xp0,xp11,xn1); 

x12span = linspace(xp11,xp12,xn2); 


xl3aspan = linspace(xp12,xcon1,xn2); 


x13bspan = linspace(xcon1,xcon2,xn2); 

xl3cspan = linspace(xcon2,xp13,xn2); 

xl4span = linspace(xp13,xp14,xn3); 

xl5span = linspace(xp14,xp15,xn2); 

xl6span = linspace(xp15,xp1,xn1); 

x2span = linspace(xp1,xp2,xn1); 

x3span = linspace(xp2,xp3,xn2); 

x4span = linspace(xp3,xpL,xn1); 

yO = yin(1:7); % make sure yO is column vector and only the first 7 elements used in a duct 
yO = yO(:); % Set up the starting boundary conditions 

[xll,yl1] = ode45(‘ductfunc2',x11span,yO); % Use the MATLAB built in ODE45 function 
yOl1 = yll(xnl,:); % Set up the boundary condition for next integration 

[x12,y12] = ode45(‘ductfunc21',x12span,y01 1); 

yO12 = yl2(xn2,:); % Set up the boundary condition for next integration 

[xl3a,yl3a] = ode45(ductfunc22al',x13aspan,y01 2); 


% From duct to a constriction, U is continuous but p,=p,-U,Z, 


gamma = 1.4; % Ratio, isobaric to isochoric specific heats 
rh = 1.2825e-2; % Hydraulic radius of the resonator, rh=r/2, m 
rhc = rh*sqrt(Rs); % Hydraulic radius in the constriction, m 


omega = 2*pi*(yl3a(xn2,6) + j.*yl3a(xn2,7));% Extract the complex frequency 
rho = 353.065./y13a(xn2,1); % Density of air as a function of temperature (K) 
c = 20.0447.*sqrt(y13a(xn2,1)); % Speed of sound in air as a function of temperature (K) 
mue = 1.846e-5.*(y13a(xn2,1)/300).%1.5.*(410.4./(y13a(xn2,1)+110.4)); % Viscosity of air 
K=2.624e2.*(y13a(xn2,1)./300).41.5.*(523.831./(y13a(xn2,1)+... 
245.4.*exp(27.6./y 13a(xn2,1)))); % Thermal conductivity, Air 
Pr = 0.60928+0.23017.*exp(-0.0028565.*y13a(xn2,1)); % Prandtl number 
dV = sqrt(2.*mue./(rho.*omega)); _% Viscous penetration depth % 
dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth % 
fv = (1-j)*dV/(2* rh); 
fve = (1-9)*dV/(2* thc); 
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fk = fv/sqrt(Pr); 

fke = fvc/sqrt(Pr); 

k = omega/c *sqrt((1+(gamma-1)*fk)/(1-fv)); | % Wave number in a regular duct 

kc = omega/c *sqrt((1+(gamma-1)*fkc)/(1-fvc)); % Wave number in the constriction 

([H,Hm0,mymz2,nynz2}=zstuff(d,h,dc,hc); 

Y l=diag(-j*sqrt(mymz2-k.%2)./(k.*rho*c)); 

Y 2=diag(-j*sqrt(nynz2-kc.%2)./(ke.*rho*c)); 

Zimag=(Hm0.'*((H* Y2*H.'+Y 1)\Hm0))/(Rs*4*1r*2); % Compute the acoustic impedance 
% caused by the constriction, using 
% higher order mode method 

Ral = rho.*omega.*dV.*(1-R)./R.*(1+(1-R.%2)./(pi*R).*log((1+R)./(1-R)))/(Rs*4*r2); 

Zal = Ral+ Zimag; 

%Zal=0; % Yf do not want to included end correction, just set the impedance to zero 

% by activating this line 

U1 = yl3a(xn2,4) + j.*y13a(xn2,5); % Volume velocity at entrance of the constriction 

yO13b = yl3a(xn2,:); 

y013b(1,2) = yO13b(1,2) - real(Zal*U1); % Correct real part of pressure 

y013b(1,3) = yO13b(1,3) - imag(Zal*U1); % Correct imag part of pressure 

[x13b,y13b] = ode45(‘ductcon’,x13bspan,y013b); % Constriction region 

% From a constriction to a duct , U is continuous but p,=p,-U,Z, 

mio = 353:065./y13b(xn2,1); 

c = 20.0447.*sqrt(y13b(xn2,1)); 

mue = 1.846e-5.*(y13b(xn2,1)/300).%1.5.*(410.4./(y13b(xn2,1)+110.4)); 

K = 2.624e-2.*(y13b(xn2,1)./300).%1.5.*(523.8306./(y13b(xn2,1)+... 

245.4.*exp(-.27.6./y13b(xn2,1))));% Thermal conductivity, Air 

Pr = 0.60928+0.23017.*exp(-0.0028565.*y13b(xn2,1));%Prandtl number 

dV = sqrt(2.*mue./(rho.*omega)); % Niscous penetration depth % 

dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth % 

fy = (1-j)*d V/(2*rh); 

fvee= (1-j)*dV/(2*rhc); 

fk = fv/sqrt(Pr); 

Pec — fve/s@rt(Pr): 

k = omega/c*sqrt((1+(gamma-1)*fk)/(1-fv)); % Wave number in regular duct 

ke = omega/c*sart((1+(gamma-1)*fkc)/(1-fvc)); % Wave number in the constriction 

Y l=diag(-j*sqrt(mymz2-k.‘2)./(k.*rho*c)); 

Y 2=diag(-j*sqrt(nynz2-ke.“2)./(ke.*rho*c)); 
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Zimag=(Hm0.'*((H* Y2*H.'+Y 1)\Hm0))/(Rs*4*r42); 

Ra2 = rho.*omega.*dV.*(1-R)./R.*(1+(1-R.%2)./(pi*R).*log((1+R)./(1-R)))/(Rs*4*r42); 
Za2 = Ra2 + Zimag; 

%La2 = O; 

U2 = y13b(xn2,4) + j*y13b(xn2,5);% Volume velocity at start of the constriction 
yOl3c = y13b(xn2,:); 

y013c(1,2) = y013c(1,2)-real(Za2*U2); % Correct real part of pressure 
y013c(1,3) = y013c(1,3)-1imag(Za2*U2); % Correct imag part of pressure 
[xl3c,y13c] = ode45(‘ductfunc22a2',x13cspan,y01 3c); 

yOISe= yl3ceGn2.)): 

[xl4,y14] = ode45(‘ductfunc23’',x14span,y013); 

yO14 = yl4(xn3,:); 

[x15,y15] = ode45(‘ductfunc24',x15span,y014); 

VOID — yldotnZ,:); 

[x16,y16] = ode45(‘ductfunc25',x16span,y0O15); 

y0O16 = yl6(xnl,:); 

[x2,y2] = ode45(‘ambhxfunc’,x2span,y016); % Integrate in the ambient heat exchanger 
% Compute the end corrections for the stack 


% Activate these command to include the end corrections in the stack region 


Yonslit = 57; % Number of the slit 

%Aslit = 2.61936e-5; % Area for one slit, m7 

%oy0=3 .06e-4; % Half spacing of the stack, m 

Yds = 4.28e-2; % Width of one slit in the stack region, m 
%hs = 6.12e-4; % Height of one slit in the stack region, m 


Jomega = 2*pi*(y2(xnl,6) + j.*y2(xnl,7)); % Complex frequency 

Jrho = 353.065./y2(xn1,1); 

Jc = 20.0447.*sqrt(y2(xnl,1)); 

Yomue = 1.846e-5.*(y2(xn1,1)/300).41.5.*(410.4./(y2(xn1,1)+110.4)); 

%dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth in regular duct 
%dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth in regular duct 
%ofv = (1-})*dV/(2* rh); 

%ofk = fv/sqrt(Pr); 

%fvs = tanh((1+j)*y0/dV)/((14j)*y0/dV); % fv for Parallel plate 

Z%ofks = fvs/sqrt(Pr); % fk for Parallel plate 

%k = omega/c *sqrt((1+(gamma-1)*fk)/(1-fv)); % Wave number in regular duct 
%oks = omega/c *sqrt((1+(gamma-1)*fks)/(1-fvs));% Wave number in the constriction 
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%(H,Hm0,mymz2,nynz2]=zstuff(d,h,ds,hs); 

% Y 1\=diag(-j*sqrt(mymz2-k.%2)./(k.*rho*c)); 

% Y 2=diag(-j*sqrt(nynz2-ks.*2)./(ks.*rho*c)); 

%Zimag=(Hm0.'*((H* Y 2*H.'+Y 1)\Hm0))/(Aslit); 

%Rast3 = rho.*omega.*dV.*(1-Rst)./Rst.*(1+(1-Rst.*2)./(pi*Rst).*log((1+Rst)./(1-... 

Rst)))/(Aslit); 
%Zastl = (Rast3 + Zimag)/nslit; % Effective impedance for end corrections of the 
% stack 

%Ustl = y2(xnl,4) + j.*y2(xnl,5); % Volume velocity at entrance of the constriction 
Zast1=0; % No end corrections for stack if this line is active 
yO2 = [y2(xnl,:) yin(8)]; % Add the 8th element (i.e. H2) in the stack region 

%yO02(1,2) = yO2(1,2) - real(Zast1*Ust1); % Correct real part of pressure 

%y02(1,3) = yO2(1,3) - imag(Zast!*Ustl); % Correct imag part of pressure 
[x3,y3] = ode45(‘stackfun',x3span,y02); % Integrate in the prime mover stack region 
% Activate these commands to include the end corrections in the stack region 

%Aslit = 2.61936e-5; % Area for one slit 

%rho = 353.065./y3(xn2,1); 

%oc = 20.0447.*sqrt(y3(xn2,1)); 

Ymue = 1.846e-5.*(y3(xn2,1)/300).%1.5.*(410.4./(y3(xn2,1)+110.4)); 

%dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth in regular duct 

%dK = sqrt(2.*K./(rho.*cp.*omega));% Thermal penetration depth in regular duct 

Jfv = (1-})*dV/(2* rh); 

%otk = fv/sqrt(Pr); 

%fvs = tanh((1+j)*yO/dV)/((1+))*yO/dV); % fv for Parallel plate 

Yotks = fvs/sqrt(Pr); % fk for Parallel plate 

%k = omega/c*sqrt((1+(gamma-1!)*fk)/(1-fv));% Wave number in regular duct 

%ks = omega/c*sqrt((1+(gamma-1)*fks)/(1-fvs));% Wave number in the constriction 

%(H,Hm0,mymz2,nynz2]=zstuff(d,h,ds,hs); 

%Y \=diag( -j*sqrt(mymz2-k.*2)./(k.*rho*c)); 

%Y 2=diag( -j*sqrt(nynz2-ks.*2)./(ks.*rho*c)); 

%Zimag=(Hm0.'*((H* Y2*H.'+Y 1)\Hm0))/(Aslit); 

%Rast4 = rho.*omega.*dV.*(1-Rst)./Rst.*(1+(1-Rst.42)./(pi*Rst).*log((1+Rst)./(1-... 

Rst)))(Aslit); | 

%Zast2 = (Rast4 + Zimag)/nslit;% Effective impedance for end corrections of stack 

Zast2=0; % If this line ic active, no end corrections for the stack is included 


%Ust2 = y3(xn2,4) + j.*y3(xn2,5);  % Volume velocity at start of the constriction 
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vO3 = y3(xn2,1:7): = % only the first 7 elements are used in the hot heat exchanger 
%y03(1.2) = yO3(1.2) - real(Zast2*Ust2); 9% Correct real part of pressure 
Gy03(1,3) = y03(1,3) - imag(Zast2*Ust2); % Correct imag part of pressure 

[x4,v4] = ode45(‘hothxfunc’,x4span,y03); 

yend = [y4(xnl.:) yin(8)]'; % This is the result at the end of the hot heat exchanger 


fvec(1.1) = yend(1,1)-Th; % Find the difference of the targeted hot temperature 
fvec(2:5,1)= yin(2:5,1)-yend(2:5,1);% Find the differences of the target pre, pim, Ure , Uim 
yend = yend(:): % Show results of the integration at the end 

% 


% End of funcvF.m 
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Function fdjacF.m: This function computes the Jacobian using finite difference. 


% Begin of function % 

function df=fdjacF(yin,fin) 

% USAGE: df=fdjacF(yin,fin) 

% 

% Inputs: y, funcvF input vector to compute Jacobian about 
Jo f evaluation of funcv at y 

% Outputs: df = matrix with df(i,j) = dfuncv(i),dy() 

% This function compute an numerical approximation to the Jacobian matrix 
% for the function funcvF 

%o 

% INPUTS: yin in the original guess vector 

Jo fin is the original differences of the target 


% Requires: funcvF.m 


% 

yin=yin(:); 

fin=fin(:); 

h=1E-10.*abs(yin)+eps; 9% make y perturbation be 1E-10 of y's current value 
% Perturb with respect to U,, U,, f,, f, and H, . i.e. y(4:8) 


% DON'T perturb with respect to Tm, pr or pi 
df=zeros(5,5); % Initialize the Jacobian matrix (5 by 5) 
for 1=4:8 

% perturb the only ith element of y % 
ypert=yin; 

ypert(i)=h(i)+yin(i); 

% find f with the perturbed y % 

pert = funcvF(ypert); 

% Forward difference formula % 

df(:,i-3) = (fpert(:)-fin(:))./h(i); % Jacobian array 
femmti(l,’."); 

end 

Jo 

% End of funcvF.m % 
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Function ductfunc2.m: This function contains the differential equations in a duct region. 


% Begin of function % 

function ydot = ductfunc2(x,y) 

% This function contains the differential equation in a dust region 

% dp,/dx = -jw*rho*U1/((1-f,)* Area) 

% dU,/dx = -jw*Area*(1+(gamma-1)*f,)*p1/(rho*c%2), U, is the volumetric velocity 
% The input y is a column vector with the component of 

%o y= ene ee Din Or Wimtletins blo): 


% in a duct region, H, 1s not used, only y(1:7) is used in duct region 
Jo 

Jo 

global T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 Th Tc r % All input temperature in C 
xpll = 0.00293; % Location of the thermocouple #9 

y(1) = Th; 

y=y(1:7); 

rh=1/2; % Hydraulic radius of the resonator 

Area = 4*r2;: % Cross-section area of the duct, m’ 

Tm = y(1); 

pri = y(2); 

pim = y(3); 

Url = y(4); 

Uim = y(5); 

lile—ay CO}. 

fim = y(7); 

U = Url + j.*Uim; % The volumetric velocity, s/m? 

p = pri + j.*pim; % Acoustic pressure 

f = frl + j*fim; % Complex frequency 

w= 27 prt: 


% Gas constants for air , Tm in Kelvin, at 1 ATM % 

% Thermal conductivity % 

K = 2.624e-2*(Tm/300)‘1.5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 
% Niscosity of Air % 

mue = 1.846e-5.*(Tm/300).%1.5.*(410.4./(Tm+110.4)); 
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gamma = 1.4; % Ratio, isobaric to isochoric specific heats 

cp=1005; % isobaric heat capacity per unit mass 

% Compute gas properties and the f function 

[fv,fk,c,rho, K, Ksolid,Prandtl]=airproperty(Tm,w, 1); 

Tmdot = (T9-Th)/xp11; % Set the temperature gradient along the duct 


% set up the ODE equation in a duct region 

temp! = -j.*w.*rho.*U./((1-fv)*Area); 

temp2 = -j.*w.*Area.*(1+(gam-1).*fk).*p./(rho.*c.42)+... 
((fk-fv)*Tmdot*U)/((1-fv)*(1-Prandtl)*Tm); 

% Set up the ordinary differential equations 

prdot = real(temp1); 

pidot = imag(temp]l); 

Urdot = real(temp2); 

Uidot = imag(temp2); 


fridot = 0; 

fimdot = 0; 

ydot =[Tmdot;prdot;pidot; Urdot; Uidot;frldot;fimdot]; 
% 


% End of ductfunc2.m % 
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Function ambhxfunc.m: This function contains the differential equations in the ambient 
heat exchanger. 


% Begin of function % 

function ydot = ambhxfunc(x,y) 

% This function contains the differential equations in the ambient heat exchanger 

% dp,/dx = -jw*rho*U1/((1-f,)* Area) 

Jo dU,/dx = -jw*Area*(1+(gamma-1)*f,)*pl/(rho*c’2), U, 1s the volumetric velocity 
% The input y is a column vector with the component of 

To va (lepe pe Ue Otek Ho); 

% in the ambient heat exchanger region, H, 1s not used 

% Note that only y(1:7) 1s used here 

% 


global T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 Th Tc % All input temperature in C 
y = y(l:7); 


y=V): % Make sure yO is a column vector 
Aambgas = 1.6547e-03; % Gas area. m’ 

yOamb = 1.7018e-03; % Half spacing of plates for ambient heat exchanger, m 
vil 1G. % Ambient heat HX at T, Kelvin 
Tm = y(1); 

pri = y(2); 

pim = y(3); 

Url = y(4); 

Uim = y(5); 

frl = y(6); 

fim = y(7); 

U = Url + j.*Uim; % The volumetric velocity, s/m° 

p = prl + j.*pim; % Acoustic pressure 

f = frl + j*fim; % Complex frequency 

Wea eo) pir ii 


% Gas properties for air , Tm in Kelvin, at 1 atm % 
% Thermal conductivity 
K = 2.624e-2*(Tm/300)41.5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 
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c = 20.0447*sqrt(Tm); % speed of sound in the air 

mue = 1.846e-5.*(Tm/300).41.5.*(410.4./(Tm+110.4)); % NViscosity of Air 
tho = 353.065/Tm; % Density of air 

gam = 1.4; % Ratio, isobaric to isochoric specific heats 
cp=1005; % Xsobaric heat capacity per unit mass 

deltaK = sqrt(2.*K./(rho.*cp.*w)); % Thermal penetration depth, m 

deltaV = sqrt(2.*mue./(rho.*w)); % Viscous penetration depth, m 
rootPrandtl=sqrt(mue.*cp./K); 

fy = tanh((1+j)*yOamb/deltaV)/((1+j)* yOamb/deltaV ); % fv for Parallel plate 
fk = fv./rootPrandtl;, % {k for Parallel plate 
% set up the ODE equations in the ambient heat exchanger region 

templ = -j.*w.*rho.*U./((1-fv)*Aambgas); 

temp2 = -j.*w.*Aambgas.*(1+(gam-1).*fk).*p./(rho.*c.%2); 


Tmdot = 0; % Note that it is assumed no temperature gradient across the ambient heat 
% exchanger 

prdot = real(temp1); 

pidot = imag(temp1); 

Urdot = real(temp2); 

Uidot = imag(temp2); 


fridot = 0; 

fimdot = 0; 

ydot =[Tmdot;prdot;pidot; Urdot; Uidot;frldot;fimdot]; 
% 


% End of ambhxfunc.m % 
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Function ductcon.m: This function contains the differential equations in the constricted 
duct. 


% Begin of function % 

function ydot = ductcon(x,y) 

% This function contains the differential equations in the constriction 

% dp,/dx = -jw*rho*U1/((1-f,)* Area) 

Jo dU,/dx = -jw*Area*(1+(gamma-1)*f,)*p1/(rho*c%2), U, is the volumetric velocity 
% The input y is a column vector with the component of 

% Y=1T pu 3PresPim3 res Vim fresfimi He] ; 

% in a duct region, H2 is not used, only y(1:7) is used in duct region 

J 


global T2 T3 T4 TS T6 T7 T8 T9 T10 T11 Th Tc Rs r 
% Define the locations of the constriction and thermoacoustics (in m) 
xpl2 = 0.084345; % Thermocouple # 10 


xp13 = 0.39615; % Thermocouple # 11 

xcon1=0.133863; % Constriction locates at xconl to xcon2 
xconz= 0:23290: % 45° long constriction, 90° from center of stack 
y=y(1:7); 

r_h = 1r/2*sqrt(Rs) ; % Hydraulic radius of the constricted duct. m 
Area = Rs*4*r42; % Cross-section area of the duct, m’ 
ey) 

prl = y(2); 

pim = y(3); 

Url = y(4); 

Uim = y(5); 

frl = y(6); 

fim = y(7); 

U = Url + j.*Uim; % Volumetric velocity, s/m° 

p = prl + j.*pim; % Acoustic pressure 

f = frl + j*fim; % Complex frequency 

w = 2*pi*f; 


% Gas properties for air , Tm in Kelvin, at 1 atm % 
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cp=1005; % isobaric heat capacity per unit mass 
c=20.0447.*sqrt(Tm); % Speed of sound in air 

nie—3 55.065.) Pin; % Density of air 

gam = 1.4; % Ratio, isobaric to isochoric specific heats 

% Thermal conductivity, Air, Tm in Kelvin 

K = 2.624e-2.*(Tm./300).%1.5.*(523.8306./(Tm+245.4.*exp(-27.6./Tm))); 
mue = 1.846e-5.*(Tm/300).%1.5.*(410.4./(Tm+110.4)); % Viscosity of Air 
Prandtl = 0.60928+0.23017.*exp(-0.0028565.*Tm); % Prandtl number 


deltaK = sqrt(2.*K./(rho.*cp.*w)); % Thermal penetration depth in the constriction, m 


deltaV = sqrt(2.*mue./(rho.*w)); % Viscous penetration depth in the constriction, m 
fee—(1-}).*deltaV /(2*r_h); % f, of the constriction 
fk = fv./sqrt(Prandtl); % f, of the constriction 


Tmdot = (T11-T10)/(xp13-xp12); % Temperature gradient across the constriction 


% set up the ODE equation in a constricted duct region 

templ = -j.*w.*rho.*U./((1-fv)*Area); 

temp2 = -j.*w.*Area.*(1+(gam-1).*fk).*p./(rho.*c.%2)+... 
((fk-fv)*Tmdot*U)/((1-fv)*(1-Prandtl)*Tm); 

prdot = real(temp1); 

pidot = imag(temp1); 


Urdot = real(temp2); 
Uidot = imag(temp2); 
fridot = 0; 
fimdot = 0; 


ydot =[Tmdot;prdot;pidot; Urdot; Uidot;frldot;fimdot]; 


% End of ductcon.m % 
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Function stackfun.m: This function contains the differential equations in the prime 


mover stack. 


% Begin of function % 

function ydot = stackfun(x,y) 

% This function contains the differential equations in the prime mover stack 
% to integrate in the stack region 

% The input y is a column vector with the component of 

To Y= lesD res Puns res Caircstinsbte!s 

Jo 

% Requires: airproperty.m 

Jo 


global T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 Th Te 


y=y(:); 

Agas = 1.7995e-3; % Total gas cross section area, m’, 56 cover glasses 
AsAgas = 0.2237; % Asolid/Agas in the stack region 

vole: % Ambient side's temperature of the prime mover stack 
Tm = y(t); 

pri = y(2); 

pim = y(3); 

Url = y(4); 

Uim = y(5); 

fitle—1y <0): 

fim = y(7); 

H2 = y(8); % Time-averaged energy flux, constant along the stack 
U = Url + j.*Uim; % Volumetric velocity, s/m’ 

p = pri + j.*pim; % Acoustic pressure 

f = frl + )*fim; % Complex frequency 

w = 2*pi*f; 


% Gas properties for air 


[fv,fk,c,rho,K,Ksolid,Prandtl] = airproperty(Tm,w,2); % index == 2 for parallel plates 


% geometry 
gam = 1.4; % Ratio, isobaric to isochoric specific heats 
cp=1005; % Isobaric heat capacity per unit mass 
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% H2 is the time-averaged energy flux, constant along the stack % 

dTmdx1 = H2/Agas-0.5/Agas*real(p*U'*(1-((fk-fv')/((1 +Prandtl)*(1-fv'))))); 

dTmdx2 = 0.5*rho*cp*abs(U)42*imag(fk+Prandtl*fv')/(real(w)*(1-Prandtl*2)*abs(1-... 
fv)*2* Agas*2) - K- AsAgas*Ksolid; 

dTmdx = dTmdx1/dTmdx2; % Temperature gradient in the stack region 

templ = -j.*w.*rho.*U./((1-fv)*Agas); 

temp2 = -j.*w.*Agas.*(1+(gam-1).*fk).*p./(rho.*c.42)+((fk-fv)*dTmdx*U)/((1-fv)*(1-... 

Prandtl)*Tm); 

Tmxdot = dTmdx; 

prdot = real(temp1); 

pidot = imag(temp1); 

Urdot = real(temp2); 

Uidot = imag(temp2); 

fridot = 0; 

fimdot = 0; 

H2dot = 0; % H2 is constant through the stack 

ydot =[Tmxdot;prdot;pidot;Urdot;Uidot;fridot;fimdot;H2dot]; 


% End of stackfun.m % 
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Function airproperty.m: This function compute the properties for air. 


% Begin of function % 
function [fv,fk,c,rho,K,Ksolid,Prandtl] = airproperty(T,w,index) 


% This function compute the air properties at 1 atm and T Kelvin 


% if index == 1 >> BL approximation for regular duct 

% index == 2 >> Parallel plate for stack region 

% index == 3 >> BL approximation for constricted duct 
global Rs 

fen. = 2.50548-2/2. % Hydraulic radius of the resonator is r/2 
yO = 3.06e-4; % Half spacing of plates for stack, m 
cp=1005; % Isobaric heat capacity per unit mass 


c=20.0447.*sqrt(T); % Speed of sound in air 
NO— bs 5.000.140 
% Thermal conductivity, Air 
K = 2.624e-2.*(T./300).%1.5.*(523.8306./(T+245.4.*exp(-27.6./T))); 
Ksolid = 0.2*(1-exp(-T/100)); % Thermal conductivity, kapton 
mue = 1.846e-5.*(T/300).41.5.*(410.4./(T+110.4)); % Viscosity of Air 
Prandtl = 0.60928+0.23017.*exp(-0.0028565.*T); % Prandtl number 
deltaK = sqrt(2.*K./(tho.*cp.*w)); % Thermal penetration depth, m 
deltaV = sqrt(2.*mue./(rho.*w)); % Viscous penetration depth, m 
if index == 
fv = (1-j).*deltaV./(2*r_h); % fv for boundary layer approximation 
elseif index == 
fv = tanh((1+j)*y0/deltaV)/((1+j)*yO0/deltaV); % fv for Parallel plate 
eiselt index == 3 
fv = (1-j).*deltaV./(2*r_h*sqrt(Rs));% fv for boundary layer approximation 
% in the constricted duct 
else 
fprintf(1,' Wrong index number !!') 
return 
end 
fk = fv./sqrt(Prandtl); 


% End of airproperty.m % 
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Function plotpUF.m: This function plot the acoustic pressure distribution in the annular 


prime mover prime mover stack 


% Begin of function % 

function [x,y,pnorm]= plotpUF(yin) 

% This function plot the acoustic pressure amplitude of the auunlar prime mover 
% It uses the solution from “updataF.m” and integrates "yin" from x=0 to x=xL 
% 


% Note: The first section of this function which is not included here, 

% contains the function ’funcvF.m”. It uses funcvF.m to integrate 

% the solution from “updateF.m” from x0 to xL 

% 

% Put :funcvF.m here 

% 

dae= Th-Tc; % Temperature difference across the prime mover stack, K 


freq = y0(6,1) +j*y0(7,1); % Complex frequency 

% Combine the solution vector 

% x 1s the position vector, y is the solution matrix 

wee sx 12:x 13a;x13b;x13e;x 143x 153x16;x2;x33x4]; 

piel (:,1:7);y 12G927)3y 13a, 177);y 130GC, Pe) sy T3c(:,1:7)3yv14C., 1:7) sy I5G, 1:7)... 
y 16(:,1:7);y2¢,1:7)sy3G, 1:7)3y4¢, 1:7); 

T= yC:,1); 

p =y(:,2)+j*y(:,3); % Complex acoustic pressure 

U =y(:,4)4j*yC,5); % Complex volume velocity 

pm=zeros(8,1); 

% Find the theoretical pressure amplitude at the locations of each microphones 


% There are 8 microphones used in the measurement 


for 1=1:8 
xm(i) = 0.03483+(i-1)*9.90375e-2; % Microphone location, m 
[Y,index] = min(abs(x-xm(i))); % Compute the corresponding index 
pc(i)=pnorm(index); % Compute the theoretical values at 
% each microphones’ locations 
end 


1S 


pc=pc(:); 


% These are the calibrated measured pressure amplitude at each microphones 


% For the constricted annular prime mover 


% Note that all these measurements are subjected to the measured temperature profile listed 


% in the function “updateF.m” 
Jo 

% Low mode, With constriction 
% Rs = 0.7, low mode 


%pm=[34.64;20.63;14.82;31.92;31.18;313.42;12.07;28.41]'; 


Yopm=[35.44;21.35153;32.63;31.94;13.78,12.16;29.07]'; 


%opm=(39.12:21.7/514.95°33.05;32:53514,0971 2.17;520089]'; 


Pomim—0.20:22.37,15.14,33.77,33.16:14.5;12.3;30.03]; 
apm =o) e276. 19052 59.00.55.242115.9/7:12.62:31.45]’; 
%pm=(37.98;324.63;15.46;36.71:36.71516.29;13.1;532.65]'; 


% Rs=0.3, low mode 
%opm=(34.89;26.09;12.32;34.67;28.18;10.54;11.6;26.91]'; 


Yopm=[39.74;29.49;14.06;38.56;31.29;11.57;13.08;29.95]'; 
%opm=[43.11;32.27515.82;42.16;34.18312.56;14.41;32.77]'; 


Yopm=[(47.39:35.68318.16;46.65;37.88;13.71516.1536.31]’; 
%pm=[58.26;45.57;24.5;60.3348.71517.44;21.35;347.11]'; 
%Yopm=(55.92;55.06;30.23;74.79;60.53;21.48;27.1358.8]'; 


% Rs=0.1, low mode 
Yopm=[40.57;37.74;34.13;38.99;28.57;10.43;9.76;26.8]'; 


Yopm=[47.45;44.12;39.87;345.48:33.24:11.91511.76;31.61]'; 


%opm=(53.45350.17345.87;52.2;38;313.43313.95;36.4]’; 
Jopm=[62.9;59.89;54.71;62.74;45.34;15.77;17.47;43.7]'; 
Ypm=[91.08;97.23;89.03;104.3;75.09;325.35;32.76;73.7]'; 
Yopm=(117.4:155.53145.9;178.6;128.7:43.3:60.38;135.1]'; 


% High mode, With constriction 

% Rs=0.7, high mode 
Jopm=[7.29;312.45;11.79;6.4734.78;9.75312.213;6.14]'; 
pm=[8.03;13.09;11.48;5.6;5.44;10.19;11.67;5.17]'; 
%pm=[8.14;12.99;11.19;5.43;5.66;10.14;11.3;4.92]’; 
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% Heater = OV 
% Heater = 8V 
% Heater = 11V 
% Heater =14V 
% Heater = 20V 
Go Heater = Za 


% Heater = OV 
% Heater =V8 
% Heater = 11V 
% Heater = 14V 
% Heater = 20V 
GanMeater—25V 


% Heater = OV 
% Heater = 8V 
% Heater = 11V 
% Heater = 14V 
Foaricater — wy 
% Heater = 25V 


% Heater = OV 
% Heater = 8V 
% Heater = 11V 


Yopm=[8.41513.22;11.19;5.39;5.97;10.3;11.23;4.84]'; % Heater = 14V 


Yopm=[8.7;13.29;10.96;5.38;6.43;10.35;310.94;4.75]'; % Heater = 20V 
epm=|8.59;13.07;10.67;5.5;6.65;10.26;10.7;4.79]'; % Heater = 25V 
% Rs=0.3, high mode 

Yopm=[6.86;12.83;6.18;6.49;2.87;9.72;11.99;4.47]'; % Heater = OV 
%Ypm=[7.24;13.41;6.43;6.56;3.05;10;12.12;4.32]'; % Heater = 8V 
%Yopm=[7.38;13.64;6.69;6.62;3.19;10.18;12.21;4.22]'; % Heater = 11V 
Yopm=[7.45;13.83;6.95;6.63;3.32;10.34;12.27;4.06]’; % Heater = 14V 
Yopm=[7.41;14.4;7.48;6.81;3.67;10.87;12.71;3.87]'; % Heater = 20V 
Yopm=[5.56;14.5;7.61;6.89;3.89;11.16;12.93;3.67]'; % Heater = 25V 
% Rs=0.1, high mode 

%pm=[8.75;15.87;14.26;8.5;2.63;11.37;12.68;4.99]'; % Heater = OV 
Yopm=(9;16.27;14.48;8.59;2.83;11.67;12.85;4.76]'; % Heater = 8V 
Ypm=[8.98;16.33314.65;8.6;2.94;11.79;13;4.56]'; % Heater = 11V 
Yopm=[8.98;16.8;15.07;8.84;3.13312.21;13.66;4.48]'; % Heater = 14) 
"aam=—[7.69;17.31;15.12;8.99;3.5:12.87;15.10;4.17]’; % Heater = 20V 
Yopm=[(6.39;17.4;14.87;9.46;3.77;13.43;15.79;4.01]'; % Heater = 25V 


% A least fit is applied to the predicted values and the measured data 
A=pm*pc/(pc'*pc); 9% Find the least square fit constant 

pmscal=pm./A; % Normalized microphone voltage by the constant A 

% These are the actual positions of the microphones, m 
xpm=[0.03483;0.13387;0.23291;0.33195;0.43099;0.53003;0.62907;0.7281 1]; 
clg 


Ix =[0.7629 0.7629]; % Mark the starting location for stack assembly 
ly = [0 1.2]; 

Ix1=[0.133863 0.133863 ]; % Mark the location of the constriction 

i= [0 1.2]; 

ez —(0.2329 0.2329 |; % Mark the location of the constriction 

i 2=(0 1.2]; 


% Plot the predicted pressure amplitude vs. the least square of the measured data 
plot(x,abs(p)./max(abs(p)),'b-',xpm,pmscal,'ro',1x,ly,'g--',1x1,ly1,'g:',1x2,ly2,'g:') 

legend(‘Matlab program ','Measured data’,’ Stack region’, 'Constriction’); 

Z%title(['Low mode with end effect 3 ; dT =',num2str(dT),’ Kelvin ;.... % Title for low mode 
Jo freq = ',num2str(freq),’ ; Rs = ',num2str(Rs)]) 

title({"'High mode with end effect 3 ; dT =',num2str(dT),’ Kelvin ;.... % Title for high mode 


Io? 


freq = ',num2str(freq),’ ; Rs = 'ynum2str(Rs)}) 
xlabel(‘Position, m') 
ylabel(’Normalized pressure’) 
axis({O 0.7923 O 1.2)) 


% End of plotpUF.m % 
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Function computeL.m: This function computes the equivalent inductance for the series 
acoustic impedance Za(w) of a size discontinuity between a duct dxh and dcxhc using 
higher order mode theory. 

This function is contributed by Ralph T. Muehleisen 


% Begin of function % 
function La=computeL(d,h,dc,hc); 


%o 

% Inputs 

% d = depth of main channel 

% h = height of main channel 

%o dc = depth of constriction 

% he = height of constriction 

% Outputs 

% La = equivalent inductance of constriction 
%o 


% This function computes the equivalent inductance for the 

% series acoustic impedance Za(w) of a size discontinuity between a 
% duct dxh and dexhc using higher order mode theory. 

% The discontinuity is assumed to be symmetric 

% in the y direction and asymmetric in the z direction 

% (i.e. going from alxbl1 to a2xb2 the open area in duct 2 is 

% (d-dc)/2<y<(d+dc)/2 and (h-hc)<z<h 

% 

% To use this add the following lines to your code 

% 

% In the beginning: 

% La = computeL(d,h,dc,hc) 

%o R=sqrt((dc*hc)./(d.*h)); 

%o Ra = rho.*w.*dV.*(1-R)./R.*( 14 (1-R.42)./(pi*R).*log((1+R)./(1-R)) ); 
% Za = Ra + j*w*La; 


R=saqrt((dc*hc)./(d.*h)); 


% compute the coupling matrix H, using n2=3 modes in small constriction 
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% and select nl based on relative sizes 
n2=s 
nl=round(n2/sqrt(R)); 


% call h2d to get the full coupling matrix 
(Hfull,my,mz,ny,nz] = h2d(R,n1,n2); 
[(MH,NH] = size(Hfull); 


HOO= Hfull (1,1); 


% extract out the higher order mode terms 
H = Hfull(2:MH,2:NH); 


% Extract out the plane wave terms 
Hm0O = Hfull(2:MH,]1); 

my = my(2:MH); 

{7 =z vila): 

ny = ny(2:NH); 

nz = nz(2:NH); 


mymz2 = (my.*pi./d).42 + (mz.*pi./h).%2; 
Mle—e pivGc), 2 + (nz.*piyne). 2; 


% The inductance approximation, it assumes k2 << my2mz2 or ny2nz2 and k~kc 
YIn = diag(sqrt(mymz2)); 

Y2n = diag(sqrt(nynz2)); 

La = rho* (Hm0.' *inv(H* Y2n*H.' + Y1in)*Hm0)/(hc*dc); 


return 


% End of computeL.m % 
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Function zstuff.m: This function computes the scattering matrices used in computation of 
series acoustic impedance Za(w) of a size discontinuity between a duct dxh and dcxhe 
using higher order mode theory. 

This function is contributed by Ralph T. Muehleisen 


% Begin of function % 

function [H,Hm0,mymz2,nynz2]=zstuff(d,h,dc,hc); 
% usage [H,Hm0,mymz?2,nynz2]=zstuff(d,h,dc,hc); 
Yo 

% Inputs 

Yo d = depth of regular duct 

Yo h = height of regular duct 


% dc = depth of constriction 

Yo he = height of constriction 

Jo 

% Outputs 

Jo H = higher order mode scattering matrix 

Yo Hm0 = scattering vector back to plane wave modes 
% mymz2, nynz2 = indices of higher order modes 
Yo 


% This function computes the scattering matrices used in computation 
% of series acoustic impedance Za(w) of a size discontinuity between a 
% duct dxh and dexhc using higher order mode theory. 

% The discontinuity is assumed to be symmetric 


%o in the y direction and asymmetric in the z direction 


Jo (i.e. going from alxbl to a2xb2 the open area in duct 2 1s 
Yo d-dc)/2<y<(d+dc)/2 and (h-hc)<z<h 
% 


% Add the following lines to your program 

Y somewhere in the beginning add 

% [H,Hm0,mymz2,nynz2)=zstuff(d,h,dc,hc); 

Jo 

% In the frequency dependent part add following. If the constriction is on the right 
%o Y l=diag( -j*sqrt( mymz2-k.42)./(k*rho*c)); 

Yo Y2=diag( -j*sqrt(nynz2-kc.*2)./(kc.*rho*c)); 
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% Zimag=(Hm0.'* ((H* Y2*H.'+Y 1)\Hm0))/Sc; 

Jo R=sqrt((dc*hc)./(d.*h)) 

Jo Ra = rho.*w.*dV.*(1-R)/R.*(14+(1-R.42)./(pi*R).*log((1+R)./(1-R))); 
Jo Za = Ra+ Zimag; 

Jo 


R=sqrt((dc*hc)./(d.*h)); 
% Compute the coupling matrix H, using n2=3 modes in small constriction 
% and select nl based on relative sizes 
n2=35; 
nl=round(n2/R); 
if nl> 10 

nl=10; 
end 
% Call h2d.m to get the full coupling matrix 
{Hfull,my,mz,ny,nz]=h2d(R,n]l ,n2); 
[MH,NH]=size(Hfull); 


HOO=Hfull(1,1); 


% Extract out the higher order mode terms 
H=Hfull(2:MH,2:NH); 


Yo extract out the plane wave terms 
Hm0 =Hfull(2:MH,]1); 
my=my(2:MH); 

mz=mz(2:MH); 

ny=ny(2:NH); 

nz=nz(2:NH); 


mymz2 = (my.*pi./d).42 + (mz.*pi./h).42; 
nynz2= (ny.*pi./dc).42 + (nz.*pi./hc).%2; 
return 


% End of zstuff.m % 
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Function h2d.m: This function is used to calculate the coupling matrix H for a 
change of size of a square duct (axa)->(bxb) with R=b/a. 
This function is contributed by Ralph T. Muehleisen 


% Begin of function % 

function [H,my,mz,ny,nz]=h2d(R,N1,N2) 
% 

% Inputs: 

% Ry = a2/al = width ratio of duct, 
% Rz = b2/b1 = height ratio of duct 


% N1= # of modes in large duct region and 

% N2 = # of modes in small duct in one direction. 
% 

% Outputs: 

% H = coupling matrix 

% Mx = large duct y mode matrix 

% Mz= large duct z mode matrix 

% Ny = small duct y mode matrix 

% Ny = small duct z mode matrix 

% 


% The routine uses the NI and N2 lowest frequency modes in each section. 
% It assumes that the y direction constriction is symmetric (i.e b is 

% centered on a) and the z constriction is asymmetric (b is at one 

% edge of a 

% 


mtZ—=N2.*2; 

Ntl=N1.42; 

% Generate ny, nz 

temp=(0:N2-1)'; 

ny=zeros(Nt2,1); 

nz=ny; 

for 1=0:N2-1 
ny(i*N2 + (temp+1),1)=[1.*ones(size(temp))]; 
nz(i*N2 + (temp+1),!)=temp; 
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end 


% now sort the modes by frequency 


f2=sqrt(ny.*2 + nz.%2); % Compute relative freq of each mode 
Ly, I]J=sort(f2); % Sort the freq. 

ny=ny(I); % Pull out ny and nz 

nz=nz(I); 


% Generate my,mz 

temp=(0:N1-1)'; 

my=zeros(Ntl,1); 

mz=my; 

for 1=0:N1-1 
my(i*N1I + (temp+1),1)=[i.*ones(size(temp))]; 
mz(i*N1 + (temp+1),1)=temp; 


end 

fl=sqrt(my.*2 + mz.‘*2); % Compute relative freq of each mode 
fy, I]=sort(f1); % Sort the freq. 

my=my(I); % Pull out ny and nz 

mz=mz(1); 


% convert my,mz,ny,nz into the matrices used to compute H=Hy*Hy; 
([Ny,My]=meshgrid(ny,my); 

[Nz,Mz]=meshgrid(nz,mz); 

Mypio2=My.*pi/2; 

Hy=zeros(size(My)); 

Hz=Hy; 

% first compute H for the symmetric y direction 

e=(1-R)/2; 

sme=sin(My*pi*e); 

cme=cos(My*pi*e); 

smepr=sin(My*pi*(e+R)); 

smep2r=sin(My*pi*(e+2*R)); 

% first fill in the general terms 
Hy=-2*R.4(3/2).*My.*(sme-(-1).“Ny.*smepr)./(pi*((My.*R).42-Ny .42+eps)); 
% find the n=m*R terms (including m=n=0) and set them to (-1)4(m+n)*Sqrt(R) 
I=find(abs(Ny-R*My)<10*eps); 
temp=(2*My.*pi.*R.*cme-sme+smep2r)./(4*My *pi*sqrt(R)+eps); 
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Hy(1)=temp(I); 

% find the m=0 terms and set them to zero 

I=find(My==0); 

Hy(1)=0*I; 

% Find the n=0 terms. they were computed above by the general expression 
% but are sqrt(2) too big 

I=find(Ny==0); 

Hy(1)=Hy(I)/sqrt(2); 

% find the m=n=0 terms and set them equal to sqrt(R) 

I=find(My==0 & Ny==0); 

Hy(I)=sqrt(R)*ones(size(I)); 

% second compute H for the asymmetric z direction 

e=(1-R); 

sme=sin(Mz*pi*e); 

cme=cos(Mz*pi*e); 

smepr=sin(Mz*pi*(e+R)); 

smep2r=sin(Mz*pi*(e+2*R)); 

% first fill in the general terms 
Hz=-2*R.4(3/2).*Mz.*(sme-(-1).4Nz.*smepr)./(p1*((Mz.*R).42-Nz.42+eps)); 
% find the n=m*R terms (including m=n=0) and set them to (-1)4(m+n)*Sqrt(R) 
I=find(abs(Nz-R*Mz)<10*eps); 
temp=(2*Mz.*pi.*R.*cme-sme+smep2r)./(4*Mz*pi*sqrt(R)+eps); 
Hz(I)=temp(1); 

% Find the m=0 terms and set them to zero 

I=find(Mz==0); 

iz t)=0*I; 

% find the n=0 terms. they were computed above by the general expression 
% but are sqrt(2) too big 

f—ind( Nz==0): 

Hz(1)=Hz(I)/sqrt(2); 

% find the m=n=0 terms and set them equal to sqrt(R) 

l=find(Mz==0 & Nz==0); 

Hz(1) = sqrt(R)*ones(size(I)); 

Be— My.*Hz; 

% End of h2d.m % 
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Function update2st.m: This function iterates the ordinary differential equation for a two- 


stack annular prime mover. 


% Begin of function % 

function [y,iterations]=update2st(y); 

Jo 

Jo 

Yo usage [newy, iterations] = update2st(y) 

Jo 

Jo 

% This function implements picards method for updating y 

% solving f(y,x)=0 for a resonator 

% Wheanoueyevector 1s.a6 Dy | vector: [1 :),.cPimi ei. if cif imiHo.3 Hoe] 
% only the last 6 elements (i.e. U,., Uin, fies fims H2,, and H,,) are updated 


% newy = y - inv(J)*f(y,x); 
% where J is the Jacobian of f(y) 
Yo REQUIRES: fdjac2st.m, func2st.m 


Nmax = 50; 
y = y(); 
iterations = 0; 
while(1) 
iterations = iterations + 1; 
fvec=func2st(y); 
checky = norm(fvec./y(1:6)) % check the norm of the ratio of the diffenence % 


if checky < 1.0E-2 % to the solution, be sure to use "/" not "/" % 
fprintf(1 ,,Done\n’); 
y=y % Print out the solution % 
[fvec,yend] = func2st(y); 
yend = yend' % print the integration results at the ends % 
Q = yend(1,6)/(2* yend(1,7)) % Find the quality factor % 
return 

elseif iterations == Nmax 


fprintf(1,"Exceed max # of iterations’); 
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end 


return 
else 
J = fdjac2st(y,fvec); % calculate the Jacobian % 
diffy = -J\fvec; % Use LU factorization to find the correction to the guess % 
if checky<1!.0e-1 % Reduce the step size when close to solution % 
daffy =diffy/2; 
elseif checky < 5.0e-2 
diffy = diffy/4; 
elseif checky < 5.0e-3 
diffy = diffy/6; 
end 
newy = y; 
newy(4:9) = y(4:9) + diffy; | % Do a Newton iteration to update y % 
y = newy % only update Ure, Uim, fre, fim, H2a, and H2b % 
fprintf(1,'Computing Update, iterations %4.0f \n',iterations) 


end 


% End of update2st.m % 
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Function fdjac2st.m: This function computes the Jacobian for the two-stack annular 


prime mover using finite difference. 


% Begin of function % 

function df=fdjac2st(yin,fin) 

% USAGE: df=fdjac2st(yin,fin) 

Jo 

% INPUTS: y funcv input vector to compute Jacobian about 

%o f evaluation of funcv at y 

% OUTPUTS: df = matrix with df(i,j) = dfuncv(i),dyq) 

% This function compute an numerical approximation to the Jacobian matrix 


% for the function funcv 


% INPUTS: 

%o yin in the original guess vector 

%o fin is the original differences of the target 
%o 


% Requires: func2st.m 


yin = yin(:); 

fin = fin(:); 

h = 1E-10.*abs(yin)+eps; % make y perturbation be 1E-10 of y's current value 
% Perturb with respect to Ur, U1, fr, fi and H2a H2b . 1.e. y(4:9) 

% DON'T perturb with respect to Tm, pr or pi 


df = zeros(6,6); % Initialize the Jacobian matrix (6 by 6) 
for I = 4:9 
% perturb the only ith element of y 

ypert=yin; 


ypert(i) = h(i)+yin(1); 
% find f with the perturbed y 
fpert = func2st(ypert); 
% Forward difference formula % 
df(:,i-3) = (fpert(:)-fin(:))./h(i); % Jacobian matrix % 
fprintf(1,".'); 
end 


% End of fdjac2st.m % 
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Function funcv2st.m: The function integrate the coupled ordinary differential equations 


of the two-stack annular prime mover. 


% Begin of function % 


function [fvec,yend] = func2st(yin) 


% This function integrates coupled ODE for the two stack annular prime move 


% Requires: 


% duct2st1.m, hothx2stl.m, stackl.m, stack2.m, ambhx2stl.m, ambhx2st2.m, 
Yo hothx2stl.m, hothx2st2.m, ode45.m 

Jo Note: Most of these functions are similar to the functions for the constricted 
% annular prime mover. Therefore, only the main functions are included here. 


global Th Tc 


% All input temperature in K 


% Before to execute the program, be sure to adjust %o 
% the target hot heat exchanger temperature to the desired value Jo 
Th = 393; 

Te = 293; 

yin = yin(:); 

xpO = 0; 

xpl = 0.16869; % Starting location for the first satck assembly 


xp2 = 0.1689948; 
xp3 = 0.1929948; 


xp4 = 0.198075; 


xp) = 0.76292; 
xp6 = 0.768; 
xp/7 = 0.792; 
maple — 0.7923; 
x= 50: 
ee 10; 

xe =—«20; 


fvec = zeros(6,1); 


% End of the first stack assembly 


% Starting location for the second satck assembly 


% End of the second stack assembly 


% Effective length of the resonator, m 


xlspan = linspace(xp0,xp1,xn1); 


x2span = linspace(xp1,xp2,xn2); 


x3span = linspace(xp2,xp3,xn3); 


x4span = linspace(xp3,xp4,xn2); 
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x5span = linspace(xp4,xp5,xn1); 

x6span = linspace(xp5,xp6,xn2); 

x7span = linspace(xp6,xp7,xn3); 

x8span = linspace(xp7,xpL,xn2); 

Jo 

% Make sure that yO is column vector, only the first 7 elements are used in the duct% 

yO = yin(1:7); 

yO = yO(:); 

[xl,yl] = ode45(‘duct2st1',x1span,y0); 

yOl = yl(xnl,:); 

y ORCI); = sib 

[x2,y2] = ode4S(‘hothx2st1',x2span,y01);% Temperature at Th 

y02 = [y2(xn2,:) yin(8)]’; % Add the 8th element (i.e. H2a) for the ist stack 

[x3,y3] = ode45(‘stack1',x3span,y02); 

¥0S = yo emioalen ) % only the first 7 elements are used in the hot heat exchanger 
[x4,y4] = ode45(‘'ambhx2st1',x4span,y03);% Target Tc at the 1st ambient heat exchanger 
y04 = y4(xn2,:); 

[x5,y5] = ode45(‘duct2st1',xSspan,y04); 

yOS = y5(xnl.,:); 


y051,1) = Te; 
[x6,y6] = ode45(‘ambhx2st2',x6span,y05);% Temperature at Tc 
y06 = [y6(xn2,:) yin(8) yin(9)]’; % Add the 8th element (i.e. H2) for the 2nd stack 


[x7,y7] = ode45(‘stack2',x7span,y06); 

yO7 = y/7(xn3,1:7); % only the first 7 elements are used in the hot heat exchanger 
[x8,y8] = ode4S(‘hothx2st2',x8span,y07); % Target Th at 2nd hot heat exchanger 

% Combine results of the integration 

% This is the result at the end of the 2nd hot heat exchanger 

yend = [y8(xn2,:) yin(8) yin(9)]'; 

% Compute the target values % 

fvec(1,1) = y03(1,1)-Tc; % Find the difference of the targeted cold temperature at ambhx2st1l 
fvec(6,1) = yend(1,1)-Th;% Find the difference of the targeted hot temperature at hothx2st2 
% Find the differences of the target pre, pim, Ure and Uim % 

fvec(2: 5.) —svin(2:5,1)-yend(2: 5,1); 

yend = yend(:); % Show results of the integration at the end 

% End of funcv2st.m % 
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APPENDIX D. PROPERTIES OF AF 45 GLASS 


Mechanical properties: 

Density at 20 °C (68 °F): p = 2.72 g/cm’ 
Young’s modulus: c = 66.0 kN/mm’ 
Torsional modulus: E = 26.7 kN/mm/ 
Poisson’s ratio: y = 0.235 

Electrical properties: 

Dielectric constant (1 MHz): €, = 6.2 


Dielectric loss factor (1 MHz): tan 6 = 9 x 10% 


Thermal properties: 


Temperature for a specific electrical resistivity of 10° Q cm: T_, = 610 °C (1130 °F) 
Viscosity and corresponding temperatures 
Temperature 


Designation 
in °C (°F) 


ear aTeD 
663 (1225) Annealing point 
883 (1621) Softening point 


Transformation temperature: T,= 662 °C (1224 °F) 


Viscosity 















log 1 (d Pas) 






Coefficient of mean linear thermal expansion: 049.39) = 4.5 X 10° K" 
Optical properties: 
Refractive indices at 20°C (68°F): n,(A = 546 nm) = 1.5275, n, (A = 546 nm) = 1.5255 


Abbe value: v, = 62.2% 


Luminous transmittance (glass thickness 1.1 mm): Typ,,., = 91.2% 
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APPENDIX E.1 CONSTRUCTION DRAWINGS FOR THE PRIME MOVER 
STACK HOLDER 


Horizontal bar 
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| 
seed 


thickness : 0.100"+0.002” 








#43 drill-through 
(for 2-56 flat head screw) 
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APPENDIX E.2 CONSTRUCTION DRAWINGS FOR THE AMBIENT 
HEAT EXCHANGER 








Vertical bar 
aaa 0.048+/-0.002” 
=. 0.047+/- 
0.200 +/- \ — |—j{ | 
0.002” __iey 0.0027 fF -0.025+/-0.002” 
ee | | 0.075+/-0.002” 
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_ > drill-threaded 
0.025+/-0.002” (for 0-80 flat head screw) 
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0.028” wide with | 0.160 “ (edge to center) 
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(center to center) | _. 0.160 “ (center to center) 
| ae 
| 
0.095 +/-0.002” 
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APPENDIX E.2(CONTINUED) CONSTRUCTION DRAWINGS FOR THE 
AMBIENT HEAT EXCHANGER 


Horizontal bar 


0.047+/-0.002” 











| — 2.030+/-0.002"— — 
a ie -: = a seer 1 —_—_—— | 
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APPENDIX E.3 CONSTRUCTION DRAWINGS FOR THE HOT HEAT 
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APPENDIX F. GRAPHS OF MATLAB AND MEASURED RESULTS 


This section which is divided into two subsections contains graphs of the results of 
the MATLAB program and measurements. First, results of the constricted annular prime 
mover are presented. Next, some predictied results for the two stack annular prime mover 


are provided. 
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APPENDIX F.1 RESULTS OF THE CONSTRICTED ANNULAR PRIME 


435 


430 


Frequency (Hz) 


410 





MOVER 


Low mode, Resonance frequency vs. DeltaT , Rs = 0.7 


— Calculated, no end corrections 
Calculated, end corrections at constriction(A) 
-—-- Calculated, end corrections at constriction(B) 
~ - Calculated, end corrections at constriction & stack (B) 
O Measured data 


50 100 150 200 
Delta T (Kelvin) 


Figure F.1-1 Comparisons of the measured and calculated resonance frequencies 


of the low mode for the constricted prime mover with the porosity of 0.7. 


Method A is the conformation transformation and Method B 1s the higher order 


mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.7 
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Figure F.1-2 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with the porosity of 0.7. 


Method A is the conformation transformation and Method B is the higher order 
mode. 
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Low mode, Resonance frequency vs. DeltaT , Rs = 0.3 
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Figure F.1-3 Comparisons of the measured and calculated resonance frequencies 
of the low mode for the constricted prime mover with the porosity of 0.3. 
Method A is the conformation transformation and Method B is the higher order 


mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.3 
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Figure F.1-4 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with the porosity of 0.3. 
Method A is the conformation transformation and Method B is the higher order 


mode. 
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Low mode with end effect 2 ; dT = 0 Kelvin ; freq = 417.9+2.149i ; Rs = 0.7 
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Figure F.1-5 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Low mode with end effect 2 ; dT = 77.03 Kelvin ; freq = 420.5+1.801i ; Rs = 0.7 
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Figure F.1-6 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=77 K. The 


calculated results are based on the higher order mode method. The frequency 1s 
the calculated value. 
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Low mode with end effect 2 ; dT = 231 Kelvin ; freq = 429.140.5389) ; Rs = 0.7 


—— Matlab program 
Measured data 


Stack region 
1} Constriction 


Normalized pressure 
oO oO 
© © 


Oo 
as 


0.2 





0 0.1 0.2 0.3 0.4 0.5 0.6 Oe, 
Position, m 


Figure F.1-7 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver 1s located 45° from the stack and AT=231 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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High mode with end effect 2 ; dT = 0 Kelvin ; freq = 438.7+5.992i : Rs = 0.7 
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Figure F.1-8 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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High mode with end effect 2 ; dT = 77.03 Kelvin ; freq = 445+6.153i ; Rs = 0.7 
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Figure F.1-9 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=77 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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High mode with end effect 2 ; dT = 231 Kelvin ; freq = 460+6.746i ; Rs = 0.7 
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Figure F.1-10 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 


calculated results are based on the higher order mode method. The frequency is 


the calculated value. 
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Low mode with end corrections ; dT = 228.2 Kelvin ; freq = 307.4-1.165i ; Rs = 0.1 
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Figure F.1-11 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=228 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. Symbol * represents the mode shape above onset, the 
measured frequency is 312 Hz. 
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Low mode with end effect 2 ; dT = 0.2 Kelvin ; freq = 361.7+2.619i : Rs = 0.3 
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Figure F.1-12 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method. The frequency is 


the calculated value. 
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Low mode with end effect 2 ; dT = 71.73 Kelvin ; freq = 362.6+1.8471 ; Rs = 0.3 


——— Matlab program 
O Measured data 
| - - Stack region 
1 : | Constriction 


Normalized pressure 
© oO 
@) 0 


Oo 


0.2 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 
Position, m 


Figure F.1-13 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=72 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Low mode with end effect 2 ; dT = 226.5 Kelvin ; freq = 369.9+0.2682i : Rs = 0.3 
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Figure F.1-14 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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High mode with end effect 2 ; dT = 0.2 Kelvin ; freq = 466.9+6.538i ; Rs = 0.3 


—— Matlab program 
O Measured data 
: -- Stack region 
1 7 Constriction 


Normalized pressure 
oO = 
o> ©o 


© 
i 


0.2 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 
Position, m 


Figure F.1-15 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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High mode with end effect 2 ; dT = 71.73 Kelvin ; freq = 468.7+6.194i : Rs = 0.3 
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Figure F.1-16 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver 1s located 45° from the stack and AT=72 K. The 


calculated results are based on the higher order mode method. The frequency is 


the calculated value. 
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High mode with end effect 2 ; dT = 226.5 Kelvin ; freq = 479.6+5.239] ; Rs = 0.3 
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Figure F.1-17 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver 1s located 45° from the stack and AT=227 K. The 


calculated results are based on the higher order mode method. The frequency 1s 
the calculated value. 
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Low mode with end effect 3 ; dT = 0 Kelvin ; freq = 295.6+3.873i ;: Rs = 0.1 
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Figure F.1-18 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Low mode with end effect 3 ; dT = 228.2 Kelvin ; freq = 306.8-0.6823i ; Rs = 0.1 
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Figure F.1-19 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=228 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 0 Kelvin ; freq = 468.8+9.638i ;: Rs = 0.1 
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Figure F.1-20 Mode shape of the high mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 228.2 Kelvin ; freq = 488.9+8.246i ; Rs = 0.1 
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Figure F.1-21 Mode shape of the high mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=228 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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Low mode with end effect 3 ; dT = 0.2 Kelvin ; freq = 361.3+2.881i : Rs = 0.3 
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Figure F.1-22 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Low mode with end effect 3 ; dT = 226.5 Kelvin ; freq = 369.6+0.4252i ; Rs = 0.3 
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Figure F.1-23 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 0.2 Kelvin ; freq = 461.9+9.625i ; Rs = 0.3 
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Figure F.1-24 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 226.5 Kelvin ; freq = 475.5+8.55i ; Rs = 0.3 
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Figure F.1-25 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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Low mode with end effect 3 ; dT = 0 Kelvin ; freq = 417.8+2.217i ;: Rs = 0.7 
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Figure F.1-26 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Low mode with end effect 3 ; dT = 231 Kelvin ; freq = 428.8+0.4826i ; Rs = 0.7 
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Figure F.1-27 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 0 Kelvin ; freq = 434.248.7791 ; Rs = 0.7 
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Figure F.1-28 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 


calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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High mode with end effect 3 ; dT = 231 Kelvin ; freq = 456.3+9.995i ; Rs = 0.7 
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Figure F.1-29 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 


calculated results are based on the higher order mode method and the end 


corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 


APPENDIX F.2 RESULTS OF THE TWO STACK ANNULAR 
PRIME MOVER 


Mode shape for of the two stack prime mover, dT = 25K; freq = 425+7.421i 
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Figure F.2-1 The calculated mode shape of the high mode of the two stack 
annular prime mover at AT = 25 K. Case 1: The duct between the two hot heat 


exchangers is held at T,. 
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Normalized pressure 


Mode shape for of the two stack prime mover, dT = -5.466e-05K; freq = 424+7.234i 
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Figure F.2-2 The calculated mode shape of the low mode of the two stack 


annular prime mover at AT = 0 K. Case 1: The duct between the two hot heat 


exchangers is held at T,. 
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